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INTRODUCTION 


Present  and  planned  Submarine  Launched  Ballistic  Missile  (SLBM)  weapons  sys¬ 
tems  employ  a  stellar  sensor  to  mitigate  some  of  the  unavoidable  error  sources 
(principally  guidance  system  alignment).  The  star  to  be  sighted  is  selected  by 
the  Fire  Control  (F/C)  program  onboard  the  submarine  and  one  of  the  criteria  used 
in  this  process  is  the  angular  proximity  of  five  of  the  planets  to  the  candidate 
stars.  If  one  of  these  brighter  planets  is  too  close  (within  some  specified  tol¬ 
erance)  to  a  particular  star,  the  observation  will  be  adversely  affected  owing  to 
the  relative  luminosities.  In  order  to  determine  the  proximity,  F/C  computes 
approximate  planetary  positions  based  on  a  simple  Keplerian  model  for  planetary 
motion  and  predetermined  sets  of  orbital  elements.  Historically  the  planetary 
elements  have  been  determined  by  least  squares  fits  to  sets  of  precise  planetary 
positions.  Because  of  supposed  accuracy  requirements,  the  fitting  spans  were  5 
years  in  duration  so  that  F/C  updates  for  orbital  elements  were  required  at  the 
same  interval. 

Recently  owing  to  an  improved  assessment  of  required  accuracies  for  plane¬ 
tary  positions,  it  has  been  recognized  that  much  longer  fit  spans  may  be  employed 
with  a  consequent  reduction  in  the  frequency  of  required  F/C  updates.  In  order  to 
facilitate  the  update  process,  a  system  of  standalone  computer  programs  has  been 
developed  to  perform  the  least  squares  fits  and  provide  quantitative  assessments 
of  the  resulting  accuracies  attained.  This  report  describes  the  formulation  which 
served  as  a  basis  for  these  programs  and  presents  the  results  of  their  initial 
exercise.  In  addition,  documentation  of  the  programs  themselves  is  contained  in 
this  report.  A  listing  of  each  program  as  well  as  sample  input  and  output  is  con¬ 
tained  in  the  appendixes.  In  addition,  Appendix  A  contains  an  outline  of  the 
program  sequence. 


ORBITAL  MODELS 


The  model  used  in  F/C  for  predicting  a  planet's  position  is  based  on  the  ap¬ 
proximation  that  the  planets  are  subject  to  only  the  central  force  gravitational 
attraction  of  the  sun.  Even  with  this  appoximation,  there  is  no  closed  form  ex- 
pression  for  planetary  positions.  Classically,  one  partitions  the  orbital  problem 
into  a  description  of  the  in-plane  evolution  of  the  planetary  position  together 
with  an  ancilliary  description  of  the  orbital  ellipse  orientation.  For  central  _  • 

force  motion,  the  planets  will  describe  elliptical  motion  with  one  focus  of  the 
ellipse  at  the  location  of  the  sun.  The  orbital  ellipse  is  completely  character¬ 
ized  by  its  semi-major  axis  (or  mean  radius)  and  its  eccentricity  (Figure  1).  The 
eccentricity  is  related  to  the  semi-major  and  semi-minor  axes  by 

e2  =  a2  -  b2  (1)  • 

a2 

1 
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FIGURE  1.  ORBITAL  ELLIPSE 


The  position  of  the  planet  in  this  orbit  is  specified  by  its  true  anomaly,  9, 
which  is  the  angular  displacement  of  the  planet  from  the  point  of  closest  approach 
(perihelion)  to  the  sun.  An  elementary  result  (Reference  1)  is  that 

r  =  aU  ,:_g2) _  (2) 

1  +  e  cos  8 
from  this  it  follows  that 

x  =  ae  +  r  cos  0  (3) 

y  =  r  sin  0  (4) 

Because  there  is  no  closed  form  expression  for  the  true  anomaly,  its  value  is 
obtained  from  a  somewhat  circuitous  procedure.  One  defines  two  intermediary  angu¬ 
lar  quantities.,  the  mean,  M,  and  eccentric,  E,  anomalies.  The  mean  anomaly  is 
just  the  average  angular  displacement  of  the  satellite  from  a  specified  epoch 

value,  M  ,  at  t  .  If  the  orbit  was  circular  (e  =  o),  then  we  would  have 
’o’  o  v  ’ 


M  =  0 


The  mean  anomaly  in  turn  is  specified  by 


M  =  M  +  n(t  -  t  ) 


NSWC  TR  83-43 


where  n,  the  mean  motion,  is  just  the  average  angular  rate  of  motion  and  it  is 

given  in  canonical  units  (see  Reference  1  and  Appendix  B  for  definition  and  con¬ 

stants)  by 

n  =  a-3^2  V 1  +  1/m  (7) 

Here  a  is  again  the  semi-major  axis  and  m  is  the  reciprocal  planet-to-sun 
mass  ratio.  (See  Appendix  C  for  planet  ratios.)  The  current  mean  anomaly,  M,  is 
used  to  obtain  the  current  eccentric  anomaly,  E,  by  an  iterative  solution  of  Kep¬ 
ler's  equation 

E  =  M  +  e  sin  E  (8) 

The  eccentric  anomaly  which  has  no  particular  interpretation  aside  from  the 
trivial  case  for  circular  orbits  can  be  used  to  obtain 

r  =  a(l  -  e  cos  E)  (9) 

r  cos  0  =  a  cos  E  -  ae  (10) 


and 


r  sin  0  =  a  "V 1  -  e2  sin  E  (11) 

thus  specifying  the  in-plane  position  of  the  planet.  In  order  to  specify  the 
planet's  position  with  respect  to  some  chosen  inertial  system,  one  needs  to  de¬ 
scribe  the  orientation  of  the  orbital  ellipse  with  respect  to  that  frame.  A 
standard  inertial  system  used  in  describing  planetary  motion  is  a  mean*  heliocen¬ 
tric,  ecliptic  frame.  This  is  a  frame  with  its  origin  at  the  sun  and  a  principle 
plane  which  approximates  the  Earth's  orbital  plane.  The  x-axis  is  specified  by 
the  intersection  of  the  ecliptic  plane  and  some  mean  equatorial  plane  with  the 
positive  direction  pointing  toward  the  vernal  equinox  which  is  in  the  general  di¬ 
rection  of  the  constellation  Aries.  Three  angular  variables  are  required  to 
specify  the  orientation  of  the  orbital  ellipse  with  respect  to  this  frame  (Figure 
2).  These  are  the  inclination,  i,  the  right  ascension  of  the  ascending  node,  ft, 
and  the  argument  of  perihelion,  U).  The  inclination  is  just  the  angle  between  the 
planet  orbital  plane  and  the  fundamental  plane.  The  intersection  of  these  two 
planes  is  called  the  line  of  nodes,  and  the  right  ascension  of  the  ascending  node 
is  just  the  angle  between  the  x-axis  and  that  part  of  the  line  of  nodes  corre¬ 
sponding  to  the  planet  ascending  through  the  ecliptic.  The  argument  of  perihelion 
is  the  in-plane  angle  between  the  line  of  nodes  and  the  direction  of  perihelion. 
Sometimes  (as  in  the  F/C  selection  program)  one  specifies  the  argument  of  latitude 
instead  of  the  argument  of  perihelion,  the  former  being  defined  as 


ui  =  ui  +  ft 


(12) 


*Various  mean  inertial  frames  corresponding  to  different  epochs  are  related  to 
each  other  by  those  transformations  describing  the  precessional  motion  of  the 
Earth's  spin  axis  (Reference  2). 
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z 


FIGURE  2.  CLASSICAL  ORBITAL  ELEMENTS 

While  these  classical  elements  have  some  appeal  because  of  their  direct  geo- 
netrical  significance  and  familiarity,  they  suffer  a  deficiency  in  that  some  of 
them  lack  definition  for  particular  orbits,  thus  Q  for  zero  inclination  orbits  or 
u  for  a  zero  eccentricity  case.  In  the  case  of  the  planets,  most  have  orbits 
which  have  rather  small  inclinations  and  most  (save  Mercury)  have  fairly  small  ec¬ 
centricities  . 


In  this  analysis,  we  choose  an  alternate  set  of  elements  to  completely  elimi¬ 
nate  the  possibility  of  any  numerical  difficulties.  Our  choice  is  the  well  known 
equinoctal  variables  defined  by  (Reference  3)  (for  posigrade  orbits). 


a  =  semi-major  axis 

h  =  e  sin  u> 

k  =  e  cos  ui 

p  =  tan(i/2)  sin  fl 

q  =  tan(i/2)  cos  Q 

\  =  uj  +  Q  +  M 
o  o 


> 


(13) 


J 


4 


NSWC  TR  83-43 


The  geometrical  significance  of  the  elements  follows:  we  first  define  the 

_  _  -> 

eccentricity  vector,  e,  as  a  vector  pointing  toward  perihelion  and  having  magni¬ 
tude  e  together  with  the  nodal  vector  N  which  points  toward  the  ascending  node  and 
has  magnitude  tan(i/2).  Then  we  define  the  equinoctal  coordinate  system  (Figure 
2)  by  choosing  the  origin  of  longitudes  as  our  x^axis  and  the  orbit  plane  normal 

as  z  .  The  origin  of  longitudes  lies  in  the  orbit  plane  and  is  displaced  from 
the  node  by  the  angle  Q.  Of  course,  yg  completes  a  right-handed  system.  With 
these  definitions,  we  provide  the  following  interpretations: 

h  =  projection  of  e  onto  y^  axis 

k  =  projection  of  e  onto  xg  axis 

p  =  projection  of  N  onto  y^  axis 

q  =  projection  of  N  onto  x^  axis 

X  =  the  mean  orbital  longitude. 

To  obtain  planetary  coordinates  given  a  set  of  equinoctal  elements,  we  pro¬ 
ceed  as  follows: 


Calculate  the  current  mean  longitude  using 


X  -  X  +  n(t  -  t  ) 
o  o' 

where  the  mean  motion  n  is  given  by 
n  =  a~3/2  Vl  +  1/ok 

where 


(canonical  units; 


(14) 


(15) 


and 


a  =  the  semi-major  axis 


m.  =  the  reciprocal  planct-to-sun  mass  ratio. 


The  eccentric  longitude  is  obtained  iteratively  from  the  Newton-Raphson  method  as 

(16) 


Fi+1 


X  +  (k  -  h F.)  sin  F.  -  (h  +  kF.)  cos  F. 

i  i  i  i 


1  -  h  sin  F.  -  k  cos  F. 

i  l 


The  orbit-plane  coordinates  are  obtained  from 
XP  =  AX  cos  F  +  B  sin  F  -  ka 
YP  =  AY  sin  F  +  B  cos  F  -  ha 


(17) 

(18) 


5 
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where 

AX  =  a(l  -  h2P) 

B  =  ahkp 

AY  =  a(l  -  k2p) 

and 


P  = 


1  +yj  1  -  h2  -  k2 

and  the  inertial  cartesian  components  are  given  by 
r  =  XPf  +  YPg 

with  the  unit  vectors  f  and  g  given  by 


f  = 


1  +  p2  +  q2 


8  = 


1  +  p2  +  q2 


1  +  q2  -  p2 


2pq 

-2p 

2pq 

1  +  p2  -  q2 


2q 


09) 

(20) 

(21) 

(22) 


(23) 


(24) 


(25) 


LEAST  SQUARES  DIFFERENTIAL  CORRECTIONS 


Now  the  true  motions  of  the  planets  are  not  governed  by  a  simple  inverse 
square  central  force  motion  so  that  our  Kepler  equations  are  only  a  crude  ap¬ 
proximation  to  the  actual  motion.  In  order  to  make  this  approximation  accurate, 
we  perform  a  least  squares  fit  for  the  Kepler  orbital  parameters  using  precise 
orbital  positions  as  "observations."  It  should  be  kept  in  mind  that  in  our  case 
the  least  squares  fit  represents  a  functionalization  process  and  not  a  minimiza¬ 
tion  of  the  effects  of  observation  noise.  Since  the  planetary  coordinates  are  not 
linearly  related  to  the  Keplerian  equinoctal  elements,  we  employ  the  least  squares 
differential  correction  method  (Reference  4) .  Thus  we  assume  that  we  have  an  ini¬ 
tial  approximation  to  our  desired  element  set  a,  h,  k,  p,  q,  A;  and  we  develop  a 
set  of  least  squares  normal  equations  for  differential  corrections  6a,  6b,  6k,  6p, 


6q,  6A.  to  this  starting  set  as  follows:  let  5  (t^)  be  our  "observations,"  tlat  is 
the  coordinates  of  the  planet  at  time  t^  obtained  from  some  precise  planetary  or¬ 
bital  model  and  let  &c(t^)  be  the  corresponding  coordinates  computed  from  our 


6 


r 


— : — -r—rr- 


i 


S 


( 
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starting  element  set  using  the  equations  presented  above.  Then  the  least 
squares  "condition  equations"  for  these  observations  are  given  by 


A.  6e  =  £  (t.)  -  3  (t.)  i=l . N 

1  1  c  i' 

where 


is  the  vector  of  diffential  corrections  based  on  our  N  observations. 


(26) 


(27) 


(A.)  = 


3R  (t.) 
c  1 

3e 


BX  (t.) 
c  1 

3a 


3X  (t.) 

i 

3h 


3X  (t.) 
c  1 

3\ 


“c<V 


8a 


8W 

3h 


aZc(tN} 

3A 


(28) 


I 


t 


is  the  matrix  of  observation  partial  derivatives.  (The  required  partials  are  pre¬ 
sented  in  Appendix  D.)  The  least  squares  normal  equations  for  the  differential 
corrections  which  minimize 

N 

£  (R(t.)  -  Rc(tt))2  (29) 

are  given  by 

ATr  A  6e  =  ATr  A  (30) 

(Tr  denotes  transpose) 


where  A  is  a  column  vector  of  "residuals"  that  is 


4  =  /  X(t.)  -  X(.(ti) 


Z(t„)  -  Zc(tjj) 

In  standard  notation,  these  equations  are  written  as 
B  6e  =  E 


(31) 


(32) 
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and  the  solution  is  obtained  as 
6e  =  B'1  E 


(33) 


Now  in  practice  one  does  uDt  actually  form  normal  equations  for  the  complete 
set  of  observations  at  once  since  the  dimensions  of  the  A  matrix  would  be  prohibi¬ 
tively  large.  Instead,  we  make  use  of  the  well-known  summability  property  of  nor¬ 
mal  equations.  In  particular,  if  we  let 

B(L)  6e  =  E  (L)  (34) 


denote  the  set  of  normal  equations  developed  for  the  subset  of  observations  (L) , 
then  the  complete  set  of  normal  equations  for  all  observations  is  given  by 


fz  B<»1 

5e  = 

'z  e«->1 

1-4 

_ 1 

i 

t-i 

i _ 

(35) 


This  allows  us  to  partition  our  "observations"  into  some  dimensionally  con¬ 
venient  subsets.  We  remark  that  because  our  observations  are  not  linearly  related 
to  our  parameter  set,  one  usually  develops  an  iterative  least  squares  solution  to 
obtain  refined  improvements  for  our  differential  corrections. 


LEAST  SQUARES  FITS 


In  the  current  analysis,  the -refined  elements  were  generated  by  a  least 
squares  fit  using  precise  planetary  positions  obtained  from  the  JPL  DE-92  ephem- 
eris  (Reference  5)  as  "observations."  Because  these  positions  were  given  in 
heliocentric  mean  equatorial  of  1950.0  coordinates,  they  were  rotated  into  the 
mean  ecliptic  system  of  1950.0  prior  to  the  fitting.  The  observations  used  were 
Cartesian  coordinates  developed  every  14  days  from  the  period  JD  2443690.5  (the 
epoch  of  the  starting  elements)  through  JD  2451544.5.  No  apriori  observations 
were  employed  and  an  implicit  standard  error  of  1  km  was  used  for  the  observa¬ 
tions.  The  starting  and  final  elements  are  given  in  Table  1.  The  semi-major 
axis,  a,  is  measured  in  km  and  the  angular  variables,  i,  to,  ft,  and  Mq,  are 
measured  in  radians. 


TABLE  1.  STARTING  AND  FINAL  ELEMENTS 


Starting  Elements 


Earth 


Mercury 


Venus 


a  149598037.4642 

e  .01671416 

i  .00007099 

u>  5.010136727 

ft  3.05571794 

M  3.07458491 

o 


57909134.07 

.20562040 

.12221041 

.50664592 

.83248134 

1  24878002 


108208902.1 

.00678807 

.05923274 

.95464052 

1.32905335 

1.16065148 
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TABLE  1.  (Cont.) 


227939292.2 

.09339831 

.03224279 

4.99540506 

.85657063 

3.89962946 


Earth _ 

149597871.9212 

.01671532949782 

.00008989551986890 

5.016964968452 

3.049977190305 

3.073477677683 


227940894.4184 

.09339472415910 

.03223176504070 

4.996203334775 

.8562002543165 

3-899288800052 


Jupiter 


778942518.0 

.04749389 

.02281176 

4.789137097 

1.74539404 

1.51508308 


Refined  Elements 


Mercut 


57909137.53653 

.2056294264117 

.1222021384653 

.5070493275543 

.8323191351502 

1.248550488749 

Jupiter 

778368220.1753 

.04825983079881 

.02280397821074 

4.795474079402 

1.746015566160 

1.505867745775 


Saturn 

1430861905.0 

.05448800 

.04342258 

5.954194377 

1.97500760 

.90983226 


_ Venus 

108208961.4686 

.006775929468209 

.05923236598127 

.9553692682276 

1.328556712384 

1.160433857370 

Saturn 

1426940733.728 

.05315110813689 

.04342907360759 

5.908659409827 

1.974220936589 

.9566213166208 


RESULTS  AND  CONCLUSIONS 


In  order  to  develop  an  assessment  of  the  accuracy  of  our  fitting  process  in 
terms  of  relevant  quantities,  we  generated  a  set  of  position  histories  for  each  of 
the  planets  of  interest  at  2-week  intervals.  These  histories  were  used  to  develop 
a  set  of  right  ascension  and  declination  values  at  the  same  interval.  These  val¬ 
ues  were  compared  with  a  similar  set  based  on  the  precise  JPL  coordinates.  The 
aggregate  statistics  are  given  in  Table  2  while  Appendix  E,  Figures  E-l  through 
E-10,  shows  the  plotted  residuals.  For  the  sake  of  comparison,  Table  3  shows 
similar  statistics  based  on  positions  developed  using  the  starting  set  of  elements 
which  are  currently  employed  in  the  selection  portion  of  the  F/C  program.  The 
plotted  residuals  corresponding  to  these  elements  are  shown  in  Figures  E-ll 
through  E-20.  It  should  be  noted  that  both  sets  of  elements,  that  is,  the  start¬ 
ing  set  and  the  improved  set  were  obtained  by  fitting  the  simple  Kepler  model  to  a 
set  of  precise  independently  generated  planetary  positions.  In  the  first  in¬ 
stance,  the  fit  span  was  5  years  in  duration,  while  in  the  latter  instance,  a  span 
of  approximately  20  years  was  employed.  The  shorter  fit  span  was  chosen  because 


9 


TABLE  2.  JPL  AND  KEPLER  COMPARISON  STATISTICS 


Residual  Error  Statistics  for  Improved  Elements  (Arc  Minutes) 


Planet 

Error 

in  Right  Ascension 

Error 

in  Declination 

Mean 

0 

Peak 

Mean 

a 

Peak 

Mercury 

.00 

.17 

.73 

.00 

.07 

-.25 

Venus 

.01 

.24 

-1.33 

-.01 

.11 

-.66 

Mars 

.09 

.53 

2.09 

-.04 

.21 

-1.15 

Jupiter 

-.01 

.64 

2.18 

-.01 

.19 

.56 

Saturn 

-.08 

1.13 

-3.85 

-.04 

.35 

1.11 

TABLE  3. 

.  CURRENTLY  USED  COMPARISON 

STATISTICS 

Residual  Error 

Statistics 

for  Starting  Elements 

(Arc  Minutes) 

Planet 

Error 

in  Right  Ascension 

Error  in  Declination 

Mean 

ct 

Peak 

Mean 

0 

Peak 

Mercury 

.29 

.14 

1.30 

-.04 

.05 

-5.8 

Venus 

.28 

.16 

1.90 

-.01 

.09 

-1.1 

Mars 

-.82 

.45 

-6.60 

-.05 

.18 

-2.96 

Jupiter 

13.7 

5.5 

45.8 

.33 

2.6 

17.5 

Saturn 

26.5 

12.5 

99.9 

7.9 

5.1 

33.4 

of  supposed  accuracy  requirements  and  because  of  the  knowledge  that  for  Saturn  and 
Jupiter  there  are  long  period  perturbations  of  the  order  of  a  few  minutes  of  arc 
in  right  ascension  and  declination.  This  of  course  is  the  reason  why  we  see 
(Appendix  E,  Figures  E-17  through  E-20)  the  secular/periodic  larger  residuals  for 
Jupiter,  Saturn,  and  to  some  extent  Mars  once  we  go  beyond  the  fit  span  used  to 
develop  the  original  set  of  elements.  The  periodic  perturbations  which  were  of 
concern  in  choosing  the  shorter  span  are  evident  (see  Appendix  E,  Figures  E-7 
through  E-10)  in  the  angular  residuals  based  on  the  20  year  fit  span.  The  latter 
set  of  orbit  parameters  are  considered  superior.  The  original  set  of  accuracy 
requirements  have  been  recognized  as  overly  stringent  which  in  turn  led  to  a 
requirement  for  a  shorter  term,  higher  accuracy  model.  Unfortunately,  this  model 
becomes  unstable  beyond  the  fit  span.  Using  the  longer  fit  span,  we  obtain  a 
representation  with  slightly  degraded  accuracy  but  which  exhibits  a  long-term 
stability  in  its  error  characteristics.  Finally,  it  is  recognized  that  if  higher 
accuracy  did  become  an  issue,  the  residual  errors  shown  could  likely  be  modeled 
with  a  small  augmentation  of  the  original  parameter  set,  principally  with  the 
addition  of  some  periodic  terms  to  the  mean  motion.  In  fact,  if  one  is  willing  to 
represent  a  large  number  of  the  observed  periodicities  in  the  principal  arguments 
of  planetary  motion,  one  can  obtain  a  very  precise  (10  seconds  of  arc  or  better) 
analytical  model  for  planetary  motion  (Reference  6) . 
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KEPLER  ORBIT  GENERATION  PROGRAM  PLANET 


The  program  PLANET  generates  planetary  positions  based  on  a  simple  Keplerian 
model  and  has  an  option  to  generate  partial  derivatives.  A  flow  chart  of  PLANET 
is  given  in  Figure  3.  The  planetary  positions  are  generated  in  the  heliocentric 
ecliptic  coordinate  system.*  The  program  PLANET  consists  of  a  driver  program 
and  subroutines  NEWT,  CONVRT,  and  PARTIAL.  The  main  program  accepts  as  an  input 
file  either  classical  or  equinoctal  elements.  If  classical  elements  are  used  as 
input,  the  equinoctal  elements  are  calculated  using  equation  (13).  If  equinoctal 
elements  are  used  as  input,  the  classical  elements  are  calculated  from  the 


inverse  transformations. 

a  =  a  (36) 

e  =  k/cos  (tan*1  (h/k))  (37) 

i  =  2  tan'1  £p/sin  (tan*1  (p/q))J  (38) 

u)  =  tan'1  (h/k)  -  tan”1  (p/q)  (39) 

Q  =  tan"1  (p/q)  (40) 

M  =  A  -  tan'1  (h/k)  (41) 

oo 


The  input  file  also  includes  an  epoch  and  ending  time  in  Julian  days  and  a  delta 
time  in  days.  The  epoch  and  end  times  are  not  necessary  for  any  calculations  and 
are  used  for  reference  only.  The  delta  time  is  the  interval  of  time  between  suc¬ 
cessive  calculations  of  the  planetary  ephemerides.  Other  input  includes  planet 
ID,  interval  between  time  lines  on  system's  output,  and  the  option  for  computing 
parital  derivatives.  The  semi-major  axis  "a"  and  delta  time  are  converted  to  can¬ 
onical  units.  The  mean  motion  is  calculated  using  equation  (7).  The  mean  motion 
is  then  used  to  calculate  current  mean  longitude  used  in  subroutine  NEWT. 

As  previously  mentioned,  the  main  program  calls  to  a  subroutine  named  NEWT. 
Subroutine  NEWT  calculates  the  eccentric  longitude  from  the  mean  longitude  at  each 
given  time  step  using  Newton's  method.  For  the  epoch  time,  the  mean  longitude,  X, 
is  equal  to  the  original  input  value.  A  starting  value  is  then  computed  for  the 
eccentric  longitude  using 

Fq  =  A.  +  k  sin  \  -  h  cos  \  +  ^ sin  2\  -  ^  cos  2\.  (42) 


*  Mean  System  of  1950.0 


to 
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FIGURE  3.  PLANET  PROGRAM  FLOW  CHART 
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and  successive  approximations  are  computed  iteratively  using  equation  (16).  The 
iteration  terminates  when  the  absolute  value  of  the  difference  between  the  current 
eccentric  longitude  and  the  previously  computed  eccentric  longitude  is  less  than 
or  equal  to  a  tolerance  (.000001745).  Control  is  then  returned  to  the  main  pro¬ 
gram. 

Main  immediately  calls  to  subroutine  C0NVRT.  This  procedure  converts  ele¬ 
ments  and  eccentric  longitude  to  inertial  Cartesian  ecliptic  positions  using  equa¬ 
tions  (17)  through  (25).  Control  again  returns  to  the  main  program  where  the 
positions  computed  in  CONVRT  are  scaled  by 


XS  =  X 

•  AU 

(43) 

YS  =  Y 

•  AU 

(44) 

ZS  =  Z 

•  AU 

(45) 

Time  is  stepped  by  adding  the  delta  time.  A  new  mean  longitude  is  computed  by 
K  =  K  +  EN  •  DT  (46) 

where 


DT  =  DELTA/ TU 


(47) 


and  EN  is  the  previously  computed  mean  motion.  A  test  is  made  to  see  if  current 
time  is  less  than  ending  time.  If  current  time  exceeds  ending  time,  the  program 
returns  to  the  beginning  to  input  a  new  set  of  data  for  another  planet.  If  cur¬ 
rent  time  is  less  than  ending  time,  the  entire  process  starting  with  subroutine 
NEWT  is  repeated  using  the  new  mean  longitude.  Termination  of  the  program  occurs 
when  the  input  value  for  planet  ID  is  greater  than  100. 

Subroutine  PARTIAL  is  a  procedure  that  may  be  invoked  by  the  correct  input 
option.  PARTIAL  computes  analytical  partial  derivatives  of  positions  with  respect 
to  the  equinoctal  elements  at  each  time  step.  Here  the  eccentric  longitude  used 
is  the  current  eccentric  longitude  at  the  current  time.  The  current  canonical 
time  from  epoch  "CT"  is  accumulated  in  main  in  the  mean  longitude  stepping  loop 


CT  =  CT  +  DT 


(48) 


The  partial  derivative  equations  used  in  this  procedure  are  contained  in  Ap¬ 
pendix  D.  The  only  quantities  relating  to  the  partial  derivative  computations 

which  change  with  time  are  X  ,  Y  ,  F,  and  CT.  Because  of  this,  many  of  the  fac- 

P  P  .. 

tors  appearing  in  the  expressions  need  only  be  computed  once.  For  example,  3f, 

3f  3g  ,  3g  ,  3q 

gp,  and  are  all  computed  once. 


Most  of  the  partial  derivatives  computed  need  to  be  scaled.  No  scaling  is 

required  for  3r. 

3a 
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The  remaining  partial  derivatives,  however,  require  scaling.  This  is  accom¬ 
plished  by  multiplying  each  one  by  a  factor  of  one  AU.  A  listing  of  the  program 
and  a  sample  of  the  input  and  output  can  be  found  in  Appendix  F. 


PROGRAM  EXTRACT 


In  order  to  perform  a  least  squares  fit,  a  high  accuracy  or  "true"  position 
history  is  needed.  This  is  provided  by  the  JPL  DE-92  ephemeris.  Because  these 
positions  are  given  in  heliocentric  mean  equatorial  of  1950.0  coordinates,  it  is 
necessary  to  rotate  them  into  the  mean  ecliptic  system  of  1950.0.  This  was  accom¬ 
plished  in  program  EXTRACT.  This  program  extracts  position  histories  from  JPL  and 
rotates  them  by  using  the  following  formula: 

E  =  .4092061859  radians  (The  mean  obliquity  of  1950.0)  (49) 


Z  /  ecliptic 


0 

cos  E 
-sin  E 


y  cos  E  +  z  sin  E 
k-y  sin  E  +  z  cos  E 


(50) 


equatorial 


(51) 


Refer  to  Appendix  G  for  a  listing  of  the  program  and  sample  input  and  output. 


LEAST  SQUARES  PROGRAM  -  LESTSQ 


To  improve  the  accuracy  of  the  preliminary  orbits,  a  set  of  improvements  to 
the  starting  elements  is  obtained  using  the  method  of  least  squares  differential 
corrections.  This  method  is  the  basis  for  program  LESTSQ.  A  flow  chart  of  LESTSQ 
is  given  in  Figure  4. 

LESTSQ  receives  its  input  from  two  files.  One  file  contains  the  high  ac¬ 
curacy  planet  coordinates  provided  by  DE-92.  The  other  file  contains  the  coor¬ 
dinates  computed  in  the  Kepler  program  as  well  as  the  partial  derivatives.  From 
this  input,  LESTSQ  computes  a  column  vector  of  residuals  using  equation  (31).  The 
residuals  are  then  used  in  conjunction  with  the  matrix  of  observation  partial 
derivatives  to  compute  the  right-hand  side  of  the  least  squares  normal  equations, 
equations  (30)  and  (32).  The  left-hand  side  of  the  least  squares  normal  equations 
or  the  normal  matrix  is  then  computed.  The  ultimate  goal  of  program  LESTSQ, 
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FIGURE  4.  LEAST  SQUARES  PROGRAM  FLOW  CHART 
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which  is  to  determine  the  vector  of  differential  corrections  based  on  the  observa¬ 
tions,  is  satisfied  by  rewriting  the  least  squares  normal  equations  into  an  equi¬ 
valent  expression  given  in  equation  (33).  A  full  listing  of  the  program  as  well 
as  sample  input  and  output  is  contained  in  Appendix  H. 


PROGRAM  GEOEQ 


The  refined  elements  obtained  using  the  least  squares  method  are  used  to 
generate  a  new  set  of  improved  position  histories  for  each  of  the  planets  of 
interest.  It  is  desired  that  these  positions  be  in  a  geocentric  equatorial  co¬ 
ordinate  frame.  Program  GEOEQ  is  written  for  this  purpose.  The  Earth's  heliocen¬ 
tric  trajectory  coordinates  are  subtracted  from  each  of  the  planet's  coordinates. 
The  planet's  coordinates  are  then  rotated  using  the  inverse  transformation  of 
equation  (51). 

The  position  histories  provided  by  JPL  are  also  rotated  using  this  program. 
A  listing  of  the  program  as  well  as  sample  input  and  output  is  contained  in  Ap¬ 
pendix  I. 


PROGRAM  DELALP 


In  order  to  determine  the  accuracy  of  the  Kepler  orbit  generation  results, 
it  is  necessary  to  compare  the  model  trajectories  with  those  of  JPL  in  terms  of 
relevant  quantities.  Program  DELALP,  see  Figure  5  for  a  flow  chart,  uses  the 
rotated  Kepler  generated  positions  and  develops  a  set  of  right  ascension  and  de¬ 
clination  values  per  trajectory  using  the  following  formulas: 

a  =  arctan  (y/x),  0  <  a  <  2  n  a  =  right  ascension  (52) 


6  =  arctan  (Z/RO) ,  RO  = 

6  =  declination 

Right  ascension  and  declination  values  are  also  generated  for  JPL  using  the  same 
formulas.  The  angles  are  returned  in  arcminutes.  This  is  accomplished  by  multi¬ 
plying  each  angle  by  a  constant  where 

constant  =  10800.00/PI  (54) 

and 

PI  =  3.1415926535898.  (55) 

The  angular  differences  are  computed  and  written  onto  an  output  file.  A  set 
of  statistics  including  maximum  and  minimum  values,  mean  values,  and  standard  de¬ 
viation  for  both  sets  of  angle  differences  is  computed  for  plotting  purposes.  The 


yjx2  +  Y2, 


n 

2 


<  6 


n 

<  2 


(53) 
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INPUT  TWO  SETS 
OF  POSITION 
HISTORIES 


IS 

THIS  END 
OF  FILE 


COMPUTE  MEAN. 
STANDARD  DEVIATION. 
MAXIMUM  8  MINIMUM 
VALUES  USING  VALUE 
IN  1 


ARE 

INITIAL  TIMES 
EQUAL 


STOP  PROGRAM 


OUTPUT  FILE 


COMPUTE  T(l)  COUNTER 
COMPUTE  RIGHT  ASCENSION 
AND  DECLINATION  FOR 
BOTH  SETS  OF  POSITION 
HISTORIES  CONVERT  INTO 
ARCMINUTES 


CALL  QUAD 
(SUBROUTINE) 

PUTS  ANGLES  IN 
CORRECT  QUADRANT 


FIGURE  5.  PROGRAM  DELALP  FLOW  CHART 
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aphed  results  are  included  in  Appendix  E,  Figures  E-l  through  E-10.  A  complete 
sting  of  the  program  as  well  as  sample  input  and  output  is  included  in  Appen- 
x  J. 


PROGRAM  SELECT 


In  order  to  complete  the  accuracy  test  of  the  Kepler  orbit  generation  re- 
ilts,  a  subroutine  was  obtained  from  a  program  called  SELECTION  which  is  the  ac- 
lal  program  used  in  F/C  on  the  submarine.  Program  SELECT  is  this  subroutine.  A 
ow  chart  of  SELECT  is  given  in  Figure  6.  SELECT  uses  data  statements  as  input 
id  computes  right  ascension  and  declination  for  all  the  planets  at  one  time.  The 
iput  elements  for  SELECT  are  slightly  different  from  the  ones  used  in  PLANET. 
1LECT  inputs  the  classical  elements  e,  i,  ft,  and  as  does  PLANET.  The  slightly 

.fferent  input  elements  are  semi-major  axis  in  astronomical  units,  argument  of 
ititude  as  defined  in  equation  (12),  and  mean  motion  as  defined  in  equation  (7). 

SELECT  steps  through  time  in  the  same  manner  as  PLANET.  It  computes  the  mean 
lomaly  for  each  planet  using  equation  (6).  It  then  calls  subroutine  ELCORD.  EL- 
)RD  converts  the  time-evolved  elements  into  Cartesian  coordinates. 

First  ELCORD  computes  an  eccentric  anomaly  iteratively.  After  computing  a 
Parting  value,  successive  approximations  to  eccentric  anomaly  are  generated.  A 
ilerance  of  .000001745  is  used  to  end  the  iteration.  The  coordinates  are  compu- 

ed  using  as  input  the  sine  and  cosine  values  for  the  angle  of  ascending  node  (ft) , 

rue  anomaly  (V),  and  the  angle  of  inclination  (i)  as  well  as  values  for  semi-ma- 

or  axis  (a) ,  eccentricity  (e) ,  the  eccentric  anomaly  (E)  and  radius  vector  (R) , 

11  at  time  "t."  The  following  are  the  equations  used  for  the  transformations. 

EO  =  Am  +  e  sin  Am  +  - — — —  -  starting  value 

c  _  r-r,  L  Am  -  E0  +  e  sin  EO  ^  , 

E  =  E0  +  - = - - eccentric  anomaly 

1  -  e  cos  EO  J 

V  =  2  arctan  (tan(E/2)  Vl  +  e,  VT^e) 

R  =  a(l  -  e  cos  E) 

X  =  R  £cos  ft  •  cos(ui  -  ft  +  V)  -  sin  ft  *  sin(w  -  ft  +  V)  •  cos  i] 

Y  =  R  £sin  ft  •  cos(ui  -  ft  +  V)  +  cos  ft  •  sin(u)  -  ft  +  V)  •  cos  ij 

Z  =  R^sin(u)  -  ft  +  V)  •  sin  ij 

The  coordinates  in  a  heliocentric  ecliptic  coordinate  system  are  then  rotated 
nto  a  heliocentric  equatorial  coordinate  system.  The  coordinates  are  transmitted 
ack  to  the  main  program  where  they  are  changed  from  heliocentric  trajectory  co- 
rdioates  to  geocentric  trajectory  coordinates.  Right  ascension  and  declination 
agles  are  then  computed  using  equations  (52)  and  (53). 


(56) 

(57) 

(58) 

(59) 

(60) 
(61 J 
(62) 
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^  START 


FIGURE  6.  PROGRAM  SELECT  FLOW  CHART 
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The  angles  are  returned  in  arc  minutes  and  written  to  five  different  output 
es  keeping  in  mind  that  right  ascension  and  declination  are  computed  for  all 
the  planets  of  interest  at  one  time.  It  is  now  necessary  to  compare  the 
ults  of  SELECT  and  JPL  using  program  DELALP.  Since  SELECT  generates  right 
ension  and  declination  angles  and  JPL  generates  position  histories,  it  is  ne- 
sary  to  change  DELALP  slightly.  An  input  statement  is  changed  to  allow  angles 
input.  The  computations  that  changed  Kepler's  position  histories  into  angles 
id  to  be  deleted.  DELALP,  renamed  DELAPM  for  clarity,  generates  a  set  of  delta 
;ht  ascension  and  delta  declinations  and  a  set  of  statistics  for  SELECT  and  JPL, 
:essary  for  plotting,  just  as  previously  done  for  Kepler  and  JPL.  The  statis¬ 
ts  are  given  in  Table  4  and  a  listing  of  DELAPM  can  be  found  in  Appendix  J. 


TABLE  4.  RESIDUAL  ERROR  STATISTICS  FOR  IMPROVED  ELEMENTS 
(ARC  MINUTES)  SELECT 


Planet  Error  in  Right  Ascension  Error  in  Declination 


Mean 

a 

Peak 

Mean 

a 

Peak 

Mercury 

.00 

.17 

.73 

.00 

.07 

-.25 

Venus 

.01 

.24 

-1.33 

-.01 

.11 

-.66 

Mars 

.09 

.53 

2.09 

-.04 

.21 

-1.15 

Jupiter 

-.01 

.64 

2.18 

-.01 

.19 

.56 

Saturn 

-.08 

1.13 

-3-85 

-.04 

.35 

1.11 

SELECT  was  implemented  to  ensure  that  the  F/C  computer,  using  the  refined 
sments  produced  by  LESTSQ,  generated  the  same  results  as  the  Kepler  orbit  gene- 
tion  program.  A  comparison  of  Table  4  and  Table  2  shows  the  two  programs  pro¬ 
ved  the  same  results.  For  this  reason,  the  graphed  results  of  program  SELECT 
s  omitted  as  they  would  be  a  duplication  of  Figures  E-l  through  E-10.  A  com- 
?te  listing  of  the  program  as  well  as  sample  input  and  output  is  given  in  Ap- 
idix  K. 


PROGRAM  PLOT 


Program  PLOT  graphs  the  differences  in  angular  positions  computed  in  Program 
tALP  as  a  function  of  time.  An  initialization  time  of  2443690.5  and  an  ending 
ne  of  2451544.5  are  used.  The  trajectories  themselves  were  computed  every  14 
ys.  Each  fifth  trajectory  point  is  plotted.  The  plotted  results  are  contained 
Appendix  E.  A  complete  listing  of  the  program  and  sample  input  is  given  in 
pendix  L,  and  a  sample  output  plot  is  shown  in  Figure  L-l. 
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APPENDIX  A 


OUTLINE  OF  PROGRAM  SEQUENCE 


Kepler  Orbit  Generation  Program  PLANET 
A.  Input 

a.  ID  for  planet 

1.  1  =  Mercury 

2.  2  =  Venus 

3.  3  =  Earth 

4.  4  -  Mars 

5.  5  =  Jupiter 

6.  6  =  Saturn 

b.  Element  option 

1.  0  =  Classical 

2.  1  =  Equinoctal 

c.  Interval  between  time  lines  on  system's  output 

d.  Option  for  computing  partial  derivatives 

1.  0  =  No  partials 

2.  1  =  Compute  partials 

e.  Epoch  time  (Julian  days) 

f.  Delta  time  (days) 

g.  End  time  (Julian  days) 

h.  Elements 

1.  Classical 

a.  a  =  Semi-major  axis 

b.  e  =  Eccentricity 

c.  i  =  Inclination 

d.  u)  =  Angle  of  perigee 

e.  =  Angle  of  ascending  node 

f.  M  =  Initial  mean  anomaly 


Equinoctal 


=  Semi-major  axis 


f.  \  -  Mean  longitude 

Output  0 

a.  Time  (Julian  days) 

b.  Position  vectors  (heliocentric  ecliptic) 

c.  Partial  derivatives  of  x,  y,  and  z  with  respect  to  equinoctal  ele¬ 
ments  . 


A- 1 
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II.  Program  EXTRACT 

A.  Input 

a.  Time  (Julian  days)  (binary) 

b.  Position  vectors  from  JPL  (heliocentric  equatorial)  (in  binary) 

B .  Output 

a.  Time  (Julian  days) 

b.  Position  vectors  (heliocentric  ecliptic) 

III.  Program  LESTSQ 

A.  Input 

a.  Time  (Julian  days) 

b.  Position  histories  (heliocentric  ecliptic)  from  Kepler  and  JPL 

c.  Observation  matrices  "H"  which  are  the  partial  derivatives  of  X,  Y, 
and  Z  with  respect  to  the  equinoctal  elements  (from  Kepler) 

B .  Output 

a.  Matrices 

Tr 

1.  H  H. 

1  l 


2.  562 


£  »Ir  "i 


i  r 

3.  BL  where  Z  is  a  matrix  of  residuals  for  each  position  his- 


4.  562 


E  "I'  2i 


b.  Parameter  improvements 

1 .  Improvement  for  a 

2.  Improvement  for  h 

3.  Improvement  for  k 

4.  Improvement  for  X 

5.  Improvement  for  q 

6.  Improvement  for  p 

Add  parameter  improvements  to  original  elements  and  rerun  PLANET. 

Program  GEOEQ 
A.  Input 

a.  Time  (Julian  days) 

b.  Position  histories  for  planets  Mercury,  Venus,  Mars,  Jupiter,  and 
Saturn  (these  histories  are  in  heliocentric  ecliptic  and  are  the 


A-2 
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ones  obtained  after  fit  (Kepler).  In  second  run,  use  JPL  position 
histories . 

c.  Position  histories  for  planet  Earth  (heliocentric  ecliptic)  from 
Kepler  after  fit  and  JPL. 

B.  Output 

a.  Time  (Julian  days) 

b.  Position  histories  for  planets  Mercury,  Venus,  Mars,  Jupiter,  and 
Saturn  in  geocentric  equatorial. 

VI .  Program  DELALP 

A.  Input 

a.  Time  (Julian  days) 

b.  Position  histories  in  geocentric  equatorial  from  Kepler  and  JPL 

B.  Output 

a.  Time  (Julian  days) 

b.  Counter 

c.  Delta  right  ascension  (arcmin)  A  =  JPL-KEPLER 

d.  Delta  declination  (arcmin)  A  =  JPL-KEPLER 

e.  Mean  for  A  right  ascension  and  A  declination 

f.  Standard  deviation  for  A  right  ascension  and  A  declination 

g.  The  maximum  values  of  A  right  ascension  and  A  declination 

h.  Minimum  values  of  A  right  ascension  and  A  declination 

VII.  Program  PLOT 

A.  Input 

a.  Title 

b.  KC  =  Number  of  characters  in  title 

c.  L  =  Column  in  which  KC  characters  will  be  written 

d.  Option  for  graphing  =  N 

1.  1  =  plot  A  right  ascension 

2.  0  =  plot  A  declination 

B .  Output 

a.  Graphs  of  A  right  ascension  and  A  declination 

/III.  Program  SELECT 

A.  Input 

a.  Data  Statements  (elements) 

1.  Semi-major  axis  in  astronomical  units 

2.  Eccentricity 

3.  Inclination 

4.  Angle  of  ascending  node 

5.  ui  longitude  of  periapses  w  +  Q 

6.  Initial  mean  anomaly 

7.  Mean  motion 

B.  Output 

a.  Angle  of  right  ascension  (arcmin) 

b.  Angle  of  declination  (arcmin) 

IX.  Change  "read"  statement  in  DELALP  to  read  in  angles  from  SELECT  instead  of 
position  histories  and  rerun  DELALP  to  obtain  A  right  ascensions  and  A  de¬ 
clinations:  A  =  JPL-SELECT. 

X.  Rerun  program  PLOT  to  obtain  graphs  of  A  right  ascensions  and  A  declina¬ 
tions  . 
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APPENDIX  B 
CANONICAL  UNITS 


Definition:  A  normalized  system  of  units  called  "canonical  units"  was 
developed  to  help  in  mathematical  calculations  since  certain  fundamental  quanti¬ 
ties  as  the  mean  distance  from  the  Earth  to  the  sun,  the  mass  and  mean  distance  of 
the  moon,  and  mass  of  the  sun  are  not  accurately  known.  It  is  assumed  that  the 
mass  of  the  sun  is  one  "mass  unit"  and  the  mean  distance  from  the  Earth  to  the  sun 
is  one  "astronomical  unit  (AU)."  We  define  the  "distance  unit  (DU)"  to  be  the 
radius  of  a  hypothetical  reference  orbit  where  the  sun  is  the  central  body.  We 
define  the  "time  unit  (TU)"  such  that  the  speed  of  a  satellite  is  1  DU/TU. 

In  PLANET 

1  AU  =  149597871.41056  km 
1  TU  =  58.13244087  day  (Reference  5) 
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APPENDIX  C 

RECIPROCAL  PLANET-TO-SUN  MASS  RATIOS 


The 

reciprocal 

M(i) 

planet-to-sun 

i=l 

Mercury 

M(l)=6023600 

i=2 

Venus 

M(2)=408523.5 

i=3 

Earth 

M(3)=332480 

i=4 

Mars 

M(4)=3098710 

i=5 

Jupiter 

M(5)=1047 . 355 

i=6 

Saturn 

M(6)=3498.5 

i=7 

Uranus 

M(7)=22869 

i=8 

Neptune 

M(8)=19314 

i=9 

Pluto 

M(9)=3000000 

C-l 
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APPENDIX  D 
PARTIAL  DERIVATIVES 


■*  3x  3v 

ar  =  Z£  f  +  JlR  S 

3a  Da  3a  8 


(D-l) 


where 


9xp  xp  3n  ,  ..  v  [hkp  cos  F  -  (l-h2P)  sin  F] 

ga  "a  2  o'  1  -  h  sin  F  -  k  cos  F 


(D-2) 


s*  -  /  -  r  ft  -  V  «> 

3x  3y 

3i  - Jf  +  _E; 

3h  3h  *  3h  8 


-  k2p)  cos  F  -  hkp  sin  F] 


(D-3) 


(D-4) 


where 


=  a  [\(fi  ♦  sia  F  -  (zhfi  +  COS  f] 

cos  F  (hkp  cos  F  -  (1  -  h2P)  sin  F) 

1  -  h  sin  F  -  k  cos  F 


CD-5) 


sin  F  +  k 


cos  F  -  1 


cos  F  [(1  -  k2p)  cos  F  -  hkP  sin  F] 
"a  1  -  h  sin  F  -  k  cos  F 


(D-6) 


*1  =  __E  f  +  _£  i 

3k  3k  3k  8 


(D— 7  ) 
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PLOTTED  RESIDUALS 


The  graphs  show  plotted  residuals  in  right  ascension  and  declination  angles 
based  on  precise  JPL  coordinates  and  coordinates  obtained  using  fitted  elements. 

Figures  E-l  through  E-10  show  plotted  residuals  using  elements  obtained 
after  an  approximate  20-year  fit  span.  Figures  E-ll  through  E-20  show  residuals 
based  on  an  old  set  of  elements  obtained  after  a  5-year  fit  span. 


MERCURY  -  RIGHT  ASCENSION 


FIGURE  E-l.  MERCURY  -  RIGHT  ASCENSION  20-YEAR  FIT  SPAN 


MERCURY  -  DECLINATION 


FIGURE  E-2.  MERCURY  -  DECLINATION  20-YEAR  FIT  SPAN 


VENUS  -  DECLINATION 
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FIGURE  E-4.  VENUS  -  DECLINATION  20-YEAR  FIT  SPAN 


JUPITER  -  RIGHT  ASCENSION 


FIGURE  E-7 .  JUPITER  -  RIGHT  ASCENSION  20-YEAR  FIT  SPAN 


JUPITER  -  DECLINATION 


FIGURE  E-8 .  JUPITER  -  DECLINATION  20-YEAR  FIT  SPAN 


SATURN  -  RIGHT  ASCENSION 
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E-10 


SATURN  -  DECLINATION 


MERCURY  -  RIGHT  ASCENSION 


FIGURE  E-ll .  MERCURY  -  RIGHT  ASCENSION  5-YEAR  FIT  SPAN 


MERCURY  -  DECLINATION 


E-13 


FIGURE  E-12.  MERCURY  -  DECLINATION  5-YEAR  FIT  SPAN 


VENUS  -  RIGHT  ASCENSION 


»  »1  "VI'.  I1’.".  I]  Iiw.  Ill  I"  1^;  y 1  '■■  '■■■  7»  '■■  •  «"  '  1  "1 
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S3inNiw  Dav 


E-14 


FIGURE  E-13.  VENUS  -  RIGHT  ASCENSION  5-YEAR  FIT  SPAN 


VENUS  -  DECLINATION 


FIGURE  E-14.  VENUS  -  DECLINATION  5-YEAR  FIT  SPAN 


MARS  -  RIGHT  ASCENSION 


(  —  I - 1 - ' r  i  i  ■  ■  0 

o  o  o  O  o  o  o  op 

i  i  i  i  •  11 


S31flNIW  DHV 


E-16 


FIGURE  E-15.  MARS  -  RIGHT  ASCENSION  5-YEAR  FIT  SPAN 


SATURN  -  RIGHT  ASCENSION 


E-20 


FIGURE  E-19.  SATURN  -  RIGHT  ASCENSION  5-YEAR  FIT  SPAN 
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PROGRAM  PLANET 

PROGRAM  PLANET( INPUT .OUTPUT . TAPES . TAPE6 ) 

PURPOSE:  TO  COMPUTE  PLANETARY  POSITIONS  AND  PARTIAL  DERIVATIVES 

BASED  ON  A  SIMPLE  KEPLER  I  AN  MODEL 
POSITIONS  AND  PART  I ALS  FOR  EACH  TIME  STEP 
INPUT:  ID=PLANET  ID 

1 “MERCURY 
2* VENUS 
3*E ARTH 
4*MARS 
5' JUPITER 
6«SATURN 
7«URANUS 
8*N£PTUNE 
9*PLUTO 
10«MOON 
1 t«SUN 

INEL:  OPTION  FOR  INPUTTING  PLANETARY  ELEMENTS 

O'CLASSICAL  ELEMENTS 
1 *EQUINOCT AL  ELEMENTS 

INC:  INTERVAL  BETWEEN  TIME  LINES  ON  SYSTEM'S  OUTPUT 

ITRAJ:  OPTION  FOR  COMPUTING  PARTIAL  DERIVATIVES 

O«=N0  PARTI  ALS 
1 *COMPUTE  PART  I ALS 
TEPOCH:  EPOCH  TIME  (JULIAN  DAY) 

CELT :  DELTA  TIME  (DAYS) 

TEND:  END  TIME  (JULIAN  DAY) 

CLASSICAL  PLANETARY  ELEMENTS  (INEL*0) 

A:  SEMI -MAJOR  AXIS  (KM) 

E:  ECCENTRICITY 

AI:  ANGLE  OF  INCLINATION  (RAD) 

COMEGA:  ANGLE  OF  ASCENDING  NODE  (RAD) 

SOMEGA:  ANGLE  OF  PERIGEE  (RAD) 

AMO:  INITIAL  MEAN  ANOMOLY  (RAD) 

EOUINOCTAL  ELEMENTS  <INEL*1) 

A:  SEMI-MAJOR  AXIS  (KM) 

H: 

AK. 

P 

0: 

ALAMO:  LATITUDE  (RAD) 

OUTPUT:  TIME,  POSITION  (X.Y.Z),  PARTIAL  DERIVATIVES  OF  X.Y.Z 

WRT  EACH  PLANETARY  ELEMENT 
SUBPROGRAMS:  NEWT.  CONVRT.  PARTIAL 

LANGUAGE :  FORTRAN 

COMMON/CONV/BETA.XP , YP , FHAT ( 3 ) , GHAT  (  3 ) 

COMMON/ELEMTS/A . H . AK , P . 0 , ALAMO 

COMMON/PART /PAR ( 18) 

DIMENSION  AM( 9 ) 

DATA  INTPE/5/, IOUT/6/ 

DATA  AU/ 14959787 1 . 4 1056/ , TU/58 . 13244087/ 

DATA  ( AM( I ) . I«1 . 9)/6023600. , 408523 . 5,332480. , 3098710. , 1047 . 355, 

1  3498.5,22869. , 19314. ,3000000./ 

REWIND  INTPE 

REWIND  I  OUT 
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90  IFIRST  =  0 

READ( INTPE .* )  ID . INEL . INC , ITRAd 
I F ( ID . GT . 100)  GO  TO  170 
READ! INTPE .• )  TEPOCH . CELT . TEND 
I F ( INEL. NEO )  GO  TO  100 
READflN.PE.*)  A.E.A1 , SOMEGA , COMEGA , AMO 
TOMEGA-SOMEGA+COMEGA 
S INO*S IN( TOMEGA ) 

COSO*COS( TOMEGA ) 

T AN1 *T  AN( AI /2  .  ) 

SO  *  SIN(COMEGA) 

CO  *  COS (COMEGA) 

H  =  E  *  S I NO 

AK*E*COSO 

P*T  ANI *  SO 

Q=TANI*CO 

ALAMO-TOMEGA+AMO 

GO  TO  1 10 

100  READ( INTPE ,* )  A , H . AK . P . 0 . ALAMO 
EMEG-ARTNQ(H, AK) 

COMEGA-ARTNQ(P.Q) 

PI *3 .  141592653898 

I  F(  EMEG.LT.  COMEGA  )EM£G*£M£G-t-2*PI 

SOMEGA* EMEG-COMEGA 

AMO- ALAMO- EMEG 

E*AK/COS( EMEG ) 

AI*2*ATAN(P/SIN(C0MEGA) ) 

1 10  CONTINUE 

DT  *  DELT/TU 
TIME  *  TEPOCH 

PRINT  2000,  ID. INEL. INC. ITRAJ, TEPOCH. DELT, TEND 
PRINT  2001 . A. E.AI .COMEGA. SOMEGA. AMO 
PRINT  2002.  A , H . AK , P , 0 , ALAMO 
A*A/AU 

EN-A  **(-1.5 ) -SORT ( 1 . ♦ 1 . /AM( ID ) ) 

ALAM* ALAMO 

CT*0 

ITER-0 

130  CALL  NEWT ( ALAM. F) 

CALL  CONVRT(F.X.Y, 2. IFIRST) 

XS*AU*X 
YS-AU-Y 
ZS-AU *Z 

WRITE  ( I OUT , 2003 )TIME.XS.YS.Z5 
IF  (ITRAd.EO.O)  GO  TO  150 

CALL  PARTIAL  ( CT . F . EN . I F IRST ) 

IFIRST  -  1 

DO  140  1*2,18 

IFU.EQ.7.0R.I.E0.13)  GOTO  140 
PAR( I )*PAR(I )*AU 
140  CONTINUE 

WRITE ( I OUT , 2008 )  PAR 
150  ITER-ITER+1 

IF  (ITER. NE. INC)  GO  TO  160 
ITER-0 

160  TIME  •  TIME+OELT 
ALAM-ALAM  +  EN-DT 
CT-CT  +  OT 

IF  (TIME. LT. TEND)  GO  TO  130 
GO  TO  90 

170  END  FILE  10UT 

C 

1000  FORMAT  (415) 

1001  FORMAT  (3E20.13) 

C 

2000  FORMAT  (-1PLANET  ID  -  -.15/*  ELEMENT  OPTION  »  *.15/ 

1  •  OUTPUT  INTERVAL  •  -.15/*  TRAJECTORY  OPT  »  *,15/ 

1  *  EPOCH  TIME  •  -.E20.13/*  DELTA  TIME  ■  -.E20.13/ 

1  •  END  TIME  •  • . E20 . 13) 

2001  FORMAT ( -OSEMI -MAJOR  AXIS  -  -.E20.13/ 

A  *  ECCENTRICITY  «  -.E20.13/ 


F-2 


n  n  n  n  n  n  n  n  non  om 


1  »  INCLINATION  =  *.E20.13/ 

2  '  •  ANG.  OF  ASCENDING  NODE  *  -.E20.13/ 

3  *  ANG.  OF  PERIGEE  *  *.£20.13/ 

4  »  INITIAL  MEAN  ANOMOLY  *  -,£20.13) 

2002  FORMAT  ( *OS£MI -MAJOR  AXIS  =  ».E20. 13/ 

1  *  H  - ,  E20.  13/*  K  =  - ,  E20.  13/*  P  *  -.E20.13/ 

A  *  0  *  * , E20. 13/ 

2  *  LATITUDE  *  *.E20.13) 

2003  FORMAT  ( 1 F 1 2 . 3 , 3E30 . 14 ) 

2004  FdRMAT(»  -.4E18.11) 

2006  FORMAT (-  TIME  X  Y-, 

1  -  Z*) 

2005  FORMAT( 1H0.3E20. 13) 

2007  FORMAT  (1H  ) 

DOS  FORMAT  (6E20.13) 

END 

SUBROUTINE  NEWT ( ALAM . F ) 

THIS  SUBROUTINE  COMPUTES  THE  ECCENTRIC  LONGITUDE  <F)  FROM  THE 
MEAN  LONGITUDE  (ALAM)  USING  NEWTON'S  METHOD 
COMMON/E  LEMTS/ A  .  H  ,  AK  ,  P  .  0 .  ALAMO 
DATA  TOL/. 1745E-5/ 

FO  *  ALAM  +  AK*SIN(ALAM)  -  H*COS(ALAM)  +  . AK-AK-H*H)/2 . -SIN 
1  ( 2 . * ALAM )  -  AK*H/2 . *COS( 2. -ALAM) 

lOO  F  -  (ALAM  ♦  (AK-H*FO)*SIN(FO)  -  (H+AK*FO)*COS( FO) )/ 

1  ( 1 . -H-SIN(FO)  -  AK»COS( FO ) ) 

IF  (ABS(F-FO).LE.TOL)  GO  TO  110 
FO*F 

GO  TO  100 
1 10  RETURN 
END 

SUBROUTINE  CONVRT  ( F , X . Y . Z . IFIRST ) 

THIS  SUBROUTING  CONVERTS  THE  EPOCH  ELEMENTS  AND  ECCENTRIC 
LONGITUDE  TO  INERTIAL  CARTESIAN  (ECLIPTIC)  ELEMENTS 

INPUT:  F.A.H.AK.P.K 

OUTPUT :  X. Y. Z.XP. YP. FHAT(3) ,GHAT( 3) 

COMMON/ELEMTS/A ,H. AK . P . 0 , ALAMO 
COMMON/CONV/BETA . XP . YP . FHAT ( 3 ) . GHAT( 3 ) 

DATA  FIRST/O/ 

IF  ( IFIRST. NE.O)  GO  TO  100 

BETA  •  1 . /( 1 . +SORT ( 1 . -H*H-AK* AK ) ) 

AX  *  A* ( 1 . -H»H*BETA) 

B  *  A*H*AK*BETA 

AY  •  A*( 1 . -AK*AK*BETA) 

CONST  *  1 ./( 1 . +P*P+Q*0) 

FHAT (  1 )  •  CONST  * ( 1 . +Q*Q-P*P ) 

FHAT ( 2 )  •  CONST  * ( 2 . *P*Q ) 

FHAT ( 3 )  *  CONST* ( -2 . *P ) 

GHAT ( 1  )  «  FHAT ( 2  ) 

GHAT ( 2 )  *  CONST *( 1 .♦P*P-0*0) 

GHAT ( 3 )  *  CONST  * ( 2 . *0) 

100  XP*AX*COS( F  )  +  B*SIN(F)  -  AK* A 
YP  •  AY *SIN( F )  ♦  B»COS( F )  -  H-A 
X  *  XP* FHAT (  1 )  ♦  YP*GHAT (  1 ) 

Y *  XP-FHAT ( 2 )  ♦  YP*GHAT(2) 

Z  «  XP*FHAT ( 3 )  ♦  YP*GHAT( 3 ) 

RETURN 

END 

SUBROUTINE  PARTIAL  ( CT . F , EN , I F IRST ) 
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I 


c 

c  this  subroutine  computes  at  each  time  step  partial 
c  derivatives  of  position  wrt  input  elements 

c 

C  INPUT:  CT  *  TIME  -  TO 

C  F  =  CURRENT  ECCENTRIC  LONGITUDE 

C  A . H , AKK . P , 0 . ALAMO , EN , B , FHAT . GHAT . XP . VP 

C  OUTPUT:  PARTIAL  DERIVATIVES  OF  POSITION  WRT  A . H , AK , P . 0 . LAMBDA 
C 

COMMON/CONV/BETA.XP. YP , FHAT (  3  ) , GHAT ( 3 ) 

COMMON/ELEMTS/  A.H.AK.P.O. ALAMO 
COMMON/PART/  PAR( 18 i 
C 

DIMENSION  PARFQ( 3 ) . PARGO( 3 ) . PARFP ( 3 ) . P ARGP ( 3 ) 

C 

IF  (IFIRST.NE.O)  GO  TO  100 
BETA  1  *  1.-BETA 
BETA3  ■  BETA--3 
AK2  *  AK* AK 
H2  *  H*H 
AK3  *  AK2* AK 
H3  *  H2*H 
P2*  P*P 
02  *  0*0 
TERM* 1 . +P2+02 
TERM2 'TERM*  *2 

PARFO(I)  *  4 . -P2*0/TERM2 
PARFQ(2)  *  2. -P*( 1 .+P2-Q2 )/TERM2 
PARFO( 3 )  *  PARFQ( 1 )/P 
PARGO( 1 )  *  PARFO( 2 ) 

PARGQ<2)  *  -PARFQ( 1 )»( 1 ,+P2)/P2 
PARG0(3)  *  PARF0(2)/P 
PARFP(I)  *  -4 . *P*( 1 .+02 )/T£RM2 
PARFP(2)  *  2 . *0* ( 1 . +02 -P2 )/TERM2 
PARFP ( 3 )  *  -2. »FHAT( 1 )/TERM 
PARGP ( 1 )  =  PARFP ( 2 ) 

PARGP (21  *  0*PARF0(1)/P 
PARGP ( 3 )  *  -PARFQ(1)/P 
100  COSF  *  COS(F) 

SINF  *  SIN(F) 

TERM1  *  1.  -  H-SINF  -  AK*COSF 

XP A  *  XP/A  -  1 . 5*EN*CT*(H»AK*BETA*C0SF  -  ( 1 . -H2*BETA ) *SINF )/TERM 1 
YPA  ■  YP/A  -  1 ,5*EN*CT*(( 1 . -AK2*BETA)*C0SF  -  H*AK*BETA*SINF )/TERM1 
XPH  *  A*(AK»(BETA  +  H2*BETA3/BETA1 )*SINF  -  ( 2 . *H*BET A  +  H3*BETA3/ 

1  BET  A  1  )  *COSF )  -  A*COSF*(H*AK*BETA*COSF  -  ( 1 . -H2*BETA )*SINF ) 

2  /TERM1 

YPH  *  A-( -AK2*H*8ETA3*SINF/BETA 1  +  AK* < 8ETA+H2 «BETA3/BETA 1 ) -COSF - 
11.)-  A-COSF* ( ( 1 . -AK2*BETA )*COSF-H*AK*BETA*SINF )/TERM1 
XPK  *  A*( -H2*AK*BETA3*C0SF/BETA1  +  H* ( BETA+AK2 *BETA3/BETA 1 ) *SINF 
1  -1.)  +  A*SINF*(H*AK*BETA*COSF  -  ( 1.-H2*BETA)*SINF) /TERM  1 
YPK  *  A-(  -(2. -AK-BETA  +  AK3*BET A3/BET A  1 ) * S I NF  +  H* ( BET  A  +AK2 

1  -BETA3/BETA 1  )*C0SF )  ♦  A • SI NF * ( (  1 . - AK2 *BET A  1 ) • COSF  -  H*AK*BETA 

2  * S I NF ) /TERM  1 

XPL  *  A*(H*AK*BETA*COSF  -  ( 1 . -H2*BETA)*SINF )/TERM1 
YPL  *  A*( ( 1 . -AK2*BETA )*COSF  -  H*AK*BETA*SINF)/TERM1 
DO  110  1*1.3 
N  *  6* ( I  -  1 )+ 1 

PAR ( N )  *  XPA  *  FHAT ( I  )  +  YPA  *GHAT ( I  ) 

PAR ( N+ 1 )  *  XPH*  FHAT ( I  )  +  YPH*GHAT ( I  ) 

PAR ( N+2 )  *  XPK  *  FHAT ( I  )  +YPK*GHAT( I  ) 

PAR ( N+3  )  *  XPL  *  FHAT ( I )  +  YPL “GHAT ( I  ) 

PARIN+4)  *  XP*PARFO(I)  ♦  YP*PARGO(I) 

PAR ( N+5 )  -  XP*PARFP(I)  +  YP*PARGP(I) 

110  CONTINUE 
RETURN 
END 
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SAMPLE  INPUT  FOR  PROGRAM  PLANET 
100=1  1  1  1 

110=2443690.500000  14.  2443788.5 

120=57909134.07000  .2001271542194  .04721092077279  .04524996816221 

130=. 041 17767074064  2.587907280000 

140=999 

100=  ID  for  planet,  element  option  (equinoctal),  time  line  interval  on  output 
partial  derivative  option  (compute  partials) 

110=  Epoch  time,  Delta  time,  end  time 

(Julian  day)  (days)  (Julian  days) 

120=  equinoctal  elements  a,  h,  k,  and  p 
130=  equinoctal  elements  q  and  A.q 

140=  next  ID  element  needed  to  terminate  program 

The  input  data  items  for  program  PLANET  are  in  free-format. 


SAMPLE  OUTPUT  FOR  PROGRAM  PLANET  FROM  PRINT  STATEMENT 


PLANET  ID  = 

ELEMENT  OPTION  = 

OUTPUT  INTERVAL  = 

TRAJECTORY  OPT  = 

EPOCH  TIME  = 

DELTA  TIME  = 

END  TIME  = 

SEMI -MAJOR  AXIS 

ECCENTRICITY 

INCLINATION 

ANG.  OF  ASCENDING  NODE 
ANG.  OF  PERIGEE 
INITIAL  MEAN  ANOMOLY 


1 

1 

1 

1 

. 2443690500000E+07 
. I400000000000E+02 
. 2443788500000E+07 

=  . 57909 13407000E+08 

=  . 2056204000000E+00 

=  . 1222104100000E+00 

=  . 8324813400000E+00 

=  . 5066459200000E+00 

=  . 1248780020000E+01 


SEMI-MAJOR  AXIS  =  . 57909 13407000E+08 

H  = 

.2001271542194E+00 

K  = 

.4721092077279E-01 

P 

. 452499681622 IE-01 

Q 

. 41 17767074064E-01 

LATITUDE 

=  . 2587907280000E+01 

SAMPLE  OUTPUT  FOR  FILE  I OUT,  PLANET  MERCURY 


CM  CO  ^  <-•  <  CM 


1=  2443690.500  - .5568851 1360619E+08  . 77232678787654E+07 

.5697 184570 1018E+07 

2=  - . 9616533255929E+00  - . 7737235343969E+07  - .6262151893769E+08  -.2030009253 
437E+08  . 6963487064908E+06  . 1071819684898E+08 

3=  . 1333687336687E+00  - . 1113041406191E+09  - . 2570707954793E+08  -.5567665816 

912E+08  - . 6330864966379E+07  -.4569143334274E+07 
4=  . 9838144986273E-01  - .8498080881027E+07  . 3563466788752E+07  -.2758438391 

076E+07  . 1538893251802E+08  . 1109616753957E+09 

5=  2443704.500  - . 50091922068032E+08  -. 43884329261 104E+08 

. 92266050 160459E+06 

6=  - . 1560402182744E+01  . 3117557983358E+08  - . 7576013329111E+08  .2684783149 

828E+08  - . 3956718373087E+07  .5439070612617E+07 

1=  Time  in  Julian  days,  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic) 
Z  (heliocentric  ecliptic) 

_  3x  3x  3x  3x  3x  3x 

3a’  3h’  3k’  3 k>  3q’  3p 

.  3y  3y  3y  3y  3y  3y 

3a’  3h’  3k’  3A.’  3q’  3p 

,  3z  3z  3z  3z  3z  3z 

3a’  3h’  3k’  3\’  3q’  3p 


5=  Time  in  Julian  days,  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic),  Z 
(heliocentric  ecliptic) 

,  3x  3x  3x  3x  3x  3x 
3a’  3h’  3k’  3X’  3q’  3p 


The  output  generated  from  the  print  statement  in  program  PLANET  has  the  following 
formats : 


1 

1 

1 


Format  (*1  Planet  ID  =  *,I5/*  Element  Option  =  *,15/ 

*0utput  interval  =  *,I5/*  Trajectory  OPT  =  *,15/ 

*Epoch  Time  =  *,E20. 13/*Delta  Time  =  *,E20.13/ 

*End  Time  =  *,E20.13) 


Format  (*0Semi-major  Axis 
A  *Eccentricity 

1  *Inclination 

*Ang.  of  Ascending  Node 
*Ang.  of  Perigee 
*Initial  Mean  Anomaly 


=*,E20. 13/ 
=*,E20. 13/ 
=*,E20 . 13/ 
=* ,E20 . 13/ 
=*,E20 . 13/ 
=*,E20. 13) 


Format  (*0Semi-major  Axis  =  *,E20.13/ 

*H  =  *,E20.13/*K  =*,E20.13/*P  =*,E20.13/ 

*Q  =  *,E20. 13/ 

*Latitude  =  *,E20.13) 
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The  output  generated  from  the  write  statement  to  file  IOUT  has  the  following 
formats.  Format  (1F12.3,  3E30.14)  for  records  one  and  five  and  Format  (6E20.13) 
for  records  two  through  four  and  record  six.  Records  one  through  four  are  repeat¬ 
ed  computing  new  trajectories  and  partials  until  ending  time  is  met.  An  end  of 
file  is  written  to  file  IOUT  terminating  the  program. 
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APPENDIX  G 


PROGRAM  EXTRACT 


PROGRAM  EXTRACT (INPUT, OUTPUT. TAPE 1.TAPE2  TAPE3 
&TAPE4.TAPE5.TAPE6.TAPE7) 

DIMENSION  EOUT (3,7), EOUTT (2,7) 

EP0CH*2443690. 5 

EPS* .4092061859 

Ml*  1 

M2*2 

M3*3 

M4*4 

M5*5 

M6*6 

continue 

READ( 1 )DA , ( ( EOUT (I. J), 1*1, 3), J*  1,7) 

I F ( EOF ( 1  ) ) 5 , 2 
CONTINUE 

IHDA.LT.  EPOCH )G0  TO  4 
DO  3  IK* 1 ,7 

EOUTT (M2. IK)*E0UT(M2. IK) 

DO  10  N* 1 , 7 

EOUT ( M2 ,N)*COS(EPS )* EOUTT (M2,N)+SIN(EPS ) *EOUT ( M3  N) 
EOUT ( M3 . N ) =COS( EPS ) -EOUT ( M3 , N ) - SIN( EPS ) * EOUTT ( M2 . N ) 
CONTINUE 

WRITE ( 2 , 40 ) ( DA , EOUT( LM.M1 ) . LM* 1 . 3  ) 

WRITE (3, 40) (DA. EOUT( LM , M2 ) , LM« 1 , 3) 

WRITE (4,40)(DA, EOUT ( LM, M3 ) , LM* 1 , 3 ) 

WRITE ( 5 . 40 )(DA,EOUT(LM, M4 ) . LM* 1 .3) 
WRITE(6.40)(DA,E0UT(LM.M5) , LM* 1 .3) 

WRITE ( 7 , 40 ) ( DA . EOUT ( LM , M6 ) . LM* 1 , 3 ) 

GO  TO  4 
CONTINUE 
ENDFILE  2 

FORMAT ( 1F12. 3, 3E30. 14) 

STOP 

END 


The  input  file  for  program  EXTRACT  is  a  binary  file.  This  is  sample  output 
for  Tape  1,  planet  "Mercury."  These  are  JPL  rotated  coordinates. 


1=  2443690.500  -.55688472186030E+08 

. 56976677220777E+07 

2=  2443704.500  -.50093107182808E+08 

. 92347938546062E+06 

3=  2443718.500  - . 11296696038796E+08 

- . 465 17635005862E+07 

4=  2443732.500  . 336586229609 16E+08 

-.7603644726349E+07 

5=  2443746.500  .53448434105633E+08 

- . 52445205866421E+07 

6=  2443760.500  . 13070481298358E+08 

. 24536453695890E+07 


. 77241282459283F+07 
- . 43883503264652E+08 
- . 68692200060998E+08 
-.54995 137 170457E+08 
- . 47048408955032E+07 
.44048817675283E+08 


G-l 


NSWC  TR  83-43 


1=  time  (Julian  days),  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic), 
Z  (heliocentric  ecliptic) 

2=  time  (Julian  days),  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic), 
Z  (heliocentric  ecliptic) 

3=  same 

4=  same 

5=  same 

6=  same 

The  output  format  for  EXTRACT  is  Format  (1F12.3,  3E30.14)  for  each  record. 
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APPENDIX  H 


PROGRAM  LESTSQ 


PROGRAM  LESTSQ! INPUT , OUTPUT , TAPES , TAPE 7 . TAPE8 ) 

TO  IMPROVE  THE  ACCURACY  OP  PRELIMINARY  ORBIT 
USING  LEAST  SOUARES  NORMAL  EOUATIONS 
DIMENSION  XS( 1000) . YS( 1000) ,ZS( 1000) ,XSS( 1000) 
DIMENSION  YSS( 1000) .ZSSl 1000) 

DIMENSION  RESIDX(3).E(6),H(3,6),BS(6,6),£S(6) 
DIMENSION  DELTEL( 6 ) , INDEX( 5 ) , TEMP ( 6 ) , B ( 6 . 6 ) 
DATA  INTPE/6/  ,  INPUT /7/ 

DATA  BS/36*0./.ES/6*0./ 

DO  100  1*1 . lOOO 

RE ADI INTPE , 250 )DA , XS ( I ) . YS( I ) , ZS( I ) 

I F ( EOF ( INTPE ) ) 1 10.5 
N*I 

RE ADI INTPE. 240) ( (H(K. J) . J*1 ,6).K*1 ,3) 

READ! INPUT. 250)DA. XSSI  I ) . YSSl I ) . ZSSl I ) 

RESIDXI 1 )*XSS( I )*XS( I ) 

RESI0X(2)*YSS( I )-YS( I ) 

RESIDX ( 3 )  *ZSS( I )'ZS(  I  ) 

DO  10  K1-1 ,6 
DO  10  IM-1 .6 
B(K1  . IM ) *0. 

DO  50  K2*  1 .6 
DO  40  IN* 1,6 
DO  30  J1*1 ,3 

B ( K2 , IN ) *B( K2 . IN )+H( U1 , K2 ) »H( J 1 , IN) 

BSIK2. IN)*BS(K2. IN)+B!K2. IN) 

CONTINUE 
WRITE! 8 , 150) 

WRITE  I  8. 240) I <B( 16. J6 ) ,06*1 .6) . 16*1 .6) 

WRITE (8.210) 

WRITE  I  8. 240) I (BSIK3, IS) , IS* 1 ,6) ,K3*1 .6) 

DO  60  IP*1 .6 
0  E I  I P ) *0 . 

DO  80  I L*  1 , 6 
DO  70  J0*1 .3 

0  El IL)*E( IL)+RESI0X(UO)*H( JO. IL ) 

0  ES(IL)*ES(IL)+E(IL) 

WRITEI8. 160) 

WRITE  I  8. 220) IE (M3) ,M3*1 ,6) 

WRITE  I  8 , 200 ) 

WRITE  I  8. 220) I  ESI  IB) . IB=1 .6) 

00  CONTINUE 
10  CONTINUE 

CALL  CROUT (0,6,0,BS,6, BB . 0 . D , INDEX . TEMP ) 
WRITE! 8 , 180) 

WRITE  I  8. 240) I (BSILM.N0) , NO* 1 .6) . LM* 1 .6) 

WRITE  I  8 . 190) 

WR I TE I  8 . 220  )  D 
DO  120  10=1.6 
20  OELTELI I0)*0. 

DO  130  M* 1 , 6 
DO  130  N* 1 ,6 

30  DELTEL(M)*DELTEL(M)+BS(M.N)*ES(N) 

WRITEI8. 170) 

WRITE! 8 . 220 IIDELTELlM) ,M* 1,6) 

50  FORMAT!*  B*HXH  TRANSPOSE**) 
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160  FORMAT ( *  £*H  TRANSPOSE  X  RESID.*-) 

170  F ORMAT ( •  THE  PARAMETER  IMPROVEMENTS  ARE-) 
180  FORMAT ( »  B  INVERSE  IS-) 

190  FORMAT ( -  THE  DETERMINANT  IS*) 

200  FORMAT ( »  THE  E  SUM  MATRIX  E S  =  - ) 

210  FORMAT ( -  THE  B  SUM  MATRIX  BS= - ) 

220  FORMAT ( E20 .  13  ) 

230  F0RMAT13E20. 13) 

240  F0RMAT(6E20. 13) 

250  FORMAT ( 1 F 1 2 . 3 , 3E30 .14) 

STOP 

END 


MPLE  INPUT  FOR  FILE  INTPE ,  PLANET  MERCURY 

This  input  was  generated  by  Kepler  planet  program. 

2443690.500  - .55688511360619E+08  . 77232678787654E+07 

.56971845701018E+07 

-.9616533255929E+00  -.7737235343969E+07  -.6262151893769E+08  -.2030009253 
437E+08  .6963487064908E+06  . 1071819684898E+08 

. 1333687336687E+00  - . 1113041406191E+0981  -.2570707954793E+08  -.5567665816 

912E+08  - . 6330864966379E+07  -.4569I43334274E+07 

.9838144986273E-01  - . 8498080881027E+07  . 3563466788752E+07  -.2758438391 

076E+07  . 153889325 1802E+08  . 1109616753957E+09 

:  2443704.500  -.50091922068032E+08  - .43884329261 104E+08 

.92266050160459E+06 

=  -. 15604021 827 44E+01  . 3117557983358E+08  - . 7576013329111E+08  .2684783149 

828E+08  - . 3956718373087E+07  .5439070612617E+07 

=  time  (Julian  days),  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic), 

Z  (heliocentric  ecliptic) 


3x 

3x 

3x 

3x 

3x 

3x 

8a’ 

3h’ 

3k’ 

3\’ 

3q* 

3p 

1 1 

dx 

§2 

§2 

§2 

3a’ 

3h’ 

3k’ 

3A* 

3q’ 

3P 

3z 

3z 

3z 

3z 

3z 

3z 

3a’ 

3h’ 

3k’ 

3A’ 

3q’ 

3p 

:  time  (Julian  days),  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic), 

Z  (heliocentric  ecliptic) 

le  input  format  for  file  INTPE  in  program  LESTSQ  is  Format  (1F12.3,  3E30.14)  for 
‘cords  one  and  five  and  Format  (6E20.13)  for  records  two  through  four  and  record 
.x. 
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LE  INPUT  FOR  FILE  INPUT,  PLANET  MERCURY 
This  input  was  obtained  from  JPL  data. 

443690.500  -.55688472186030E+08 
.56976677220777E+07 

1443704.500  - .50093107182808E+08 
. 92347838546062E+06 

1443718.500  - . 11296696038796E+08 
- . 465 1 7635005862E+07 

’443732.500  . 336586229609 16E+08 

- . 76036447262349E+07 
>443746.500  .53448434105633E+08 

- . 52445205866421E+07 
>443760 . 500  . 13070481298358E+08 

. 24536453695890E+07 

time  (Julian  Day),  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic), 
l  (heliocentric  ecliptic) 

time  (Julian  Day),  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic), 
Z  (heliocentric  ecliptic) 

time  (Julian  Day),  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic), 
Z  (heliocentric  ecliptic) 

same 

same 


. 77241282459283E+07 
- . 43883503264652E+08 
- . 68692200060998E+08 
-.54995137170457E+08 
- . 47048408955032E+07 
.44048817675283E+08 


same 


input  format  for  file  INPUT  in  program  LESTSQ  is  Format  (1F12.3,  3E30.14)  for 
h  record. 


PLE  OUTPUT  FOR  PROGRAM  LESTSQ 


B=HXH  TRANSPOSE^ 

.9522432474214E+00  - . 8240007705156E+67  . 5714225032370E+08 

33E+08  - . 1713633537292E-06  . 2801418304443E-05 

-. 8240007705 156E+07  . 1252069390838E+17  . 3315539197539E+16 

53E+16  . 5684872774240E+15  - .5173259315036E+15 

. 5714225032370E+08  . 3315539197539E+16  . 4595006868502E+16 

99E+16  . 1739795855295E+15  - . 1583221906994E+15 

. 1182474693233E+08  .6377550616053E+16  . 2692677305999E+16 

42E+16  . 2958960392146E+15  -.2692667016374E+15 

- . 1713633537292E-06  .5684872774240E+15  . 1739795855295E+15 

46E+15  . 2773839967877E+15  . 1743971966724E+16 


.1182474693 

.6377550616 

.2692677305 

.3519593004 

.2958960392 
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Tr 

E  =  Hir  •  Z 

100=  . 14481 19382968E+05 

101=  .2751285514654E+12 

102=  . 2601902869096E+12 

103=  - . 1730909741358E+12 
104=  - . 2609827959312E+1 1 
105=  . 1312392486448E+11 

106=  THE  E  SUM  MATRIX  ES= 
107=  . 2318293753911E+05 

108=  .5366959346823E+11 

109=  . 6037889790544E+12 

110=  - . 4007279949659E+12 


300=  829761039 1578E+10 

301=  .4121374502542E+11 

302=  THE  E  SUM  MATRIX  ES= 

303=  . 2167388785937E+06 

304=  . 1006521905542E+13 

305=  .3678583926613E+12 

306=  - . 1540834753985E+13 

307=  - . 3845430331050E+12 
308=  .4075121829653E+12 

309=  B=HXH  TRANSPOSE= 

310=  . 4028520717956E+03  . 137572828796E+10  -.2033913069406E+10 

994E+10  - . 8430494574897E+08  . 7671787271573E+08 


-.1407023890 


The  output  format  for  the  normal  matrix  and  the  sum  of  the  normal  matrices  is 
Format  (6E20.13)  as  in  records  two  through  six.  The  format  for  the  right-hand 
side  of  the  normal  equations  is  Format  (E20.13)  as  in  records  100  through  105 


1000= 

1001= 

1002= 

1003= 

1004= 

1005= 

1006= 

1007= 

1008= 

1009= 

1010= 


E  =  HIr 


- . 5438037886500E+10 
.2185069505625E+11 
THE  E  SUM  MATRIX  ES= 
. 1071446905724E+07 
.  1972181680237E+13 
. 1028334456689E+12 
- .  2978376667 17 1E+13 
- .  9268345093082E+12 
.9422572290006E+12 
B=HXH  TRANSPOSE= 

. 4156773953167E+04 


.5442250291893E+10  -.5665462231553E+10  -.4450849843 


865E+10  - . 2759070951785E+09  .2510766743425E+09 


15746=  THE  PARAMETER  IMPROVEMENTS  ARE 

15747=  . 3466525731675E+01 

15748=  .2016734654506E-04 

15749=  -.4620221768938E-04 

15750=  .ll6',i45388865E-04 

15751=  . 4544758284259E-05 

15752=  -. 97496376358 15E-05 

15753: {EOR  } 
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15746=  THE  PARAMETER  IMPROVEMENTS  ARE 
15747=  6  improvement  for  "a" 

15748=  6  improvement  for  "h" 

15749=  6  improvement  for  "k" 

15750=  6  improvement  for  "A" 

15751=  6  improvement  for  "q" 

15752=  6  improvement  for  "p" 

The  output  format  for  the  parameter  improvements  is  Format  (E20.13)  as  in  records 
15746  through  15752. 
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E 

APPENDIX  I 
PROGRAM  GEOEQ 


PROGRAM  GEOEQ! INPUT . OUTPUT . TAPE6 , TAPE9 , TAPE  10) 
EPS* . 409206 1859 
4  CONTINUE 

READ16. 1 JDA.X.V ,2 
I F ( EOF ( 6 ) )40 , 10 
10  READC9. 1 IDA.XS.VS.ZS 

XX*X-XS 
YY*Y - YS 
ZZ=Z-ZS 

YN=COS( EPS  )*YY-SIN( EPS)*ZZ 
ZN*SIN( EPS  )*YY+COS( EPS )*Z2 
WRITE! 10. 1 JDA.XX.YN.ZN 
GO  TO  4 

40  CONTINUE 

ENDFILE  10 

1  FORMAT! 1F12.3.3E30. 14) 

STOP 

END 


SAMPLE  INPUT  FOR  PROGRAM  GEOEQ 

This  is  input  for  Tape  6  planet  "Mercury."  Tape  6  will  change  according  to 
the  planet.  The  planets  used  are:  Mercury,  Venus,  Mars,  Jupiter,  and  Saturn. 


1=  2443690.500  - .556861 15732592E+08 

.56958048221782E+07 

2=  2443704.500  -.50087555986015E+08 

.92061715596236E+06 

3=  2443718.500  - . 11288687930139E+08 

-.46539081624056E+07 

4=  2443732.500  . 33664975344491E+08 

-. 760393886821 13E+07 

5=  2443746.500  .53450776728722E+08 

- . 52432159266678E+07 

6=  2443760.500  . 13071876612774E+08 

. 24544754679073E+07 


. 77215772477416E+07 
- . 43887541095688E+08 
- . 68693434188216E+08 
-.54992892344885E+08 
-.47027915944105E+07 
. 44048605506 134E+08 


1=  time  (Julian  day),  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic), 
Z  (heliocentric  ecliptic) 
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2=  time  (Julian  day),  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic), 
Z  (heliocentric  ecliptic) 

3=  same 

4=  same 


5=  same 
6=  same 


The  input  format  for  Tape  6,  GEOEQ  is  Format  (1F12.3,  3E30.14)  for  each  record. 


Input  Tape  9  is  reserved  for  planet  Earth  only. 


SAMPLE  INPUT  FOR  TAPE  9  PLANET  EARTH 


1=  2443690.500  . 222884092905 18E+08 

. 13284830829530E+05 

2=  2443704.500  .56409568937832E+08 

.12177197900465E+05 

3=  2443718.500  .87428563942007E+08 

. 1040005 1418654E+05 

4=  2443732.500  . 11362338791314E+09 

. 80492806819268E+04 

5=  2443746.500  . 1335 121 1082354E+09 

. 52520508739849E+04 

6=  2443760.500  . 14593133469083E+09 

. 2161294891325 7E+04 


-. 15045085 103602E+09 
- . 14121242573654E+09 
- . 12420994559592E+09 
- . 10035637247911E+09 
- . 70936103646182E+08 
- . 37550648582708E+08 


1=  time  (Julian  days),  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic), 
Z  (heliocentric  ecliptic) 


2=  time  (Julian  days),  X  (heliocentric  ecliptic),  Y  (heliocentric  ecliptic), 
Z  (heliocentric  ecliptic) 

3-6=  same 


The  input  format  for  planet  Earth’s  GEOEQ  is  Format  (1F12.3,  3E30.14)  for  each 
record. 


SAMPLE  OUTPUT  FOR  TAPE  10 

The  following  are  the  rotated  Kepler  coordinates  for  planet  Mercury. 

1=  2443690.500  -. 779745250231 10E+08  . 14285226329953E+09 

. 68147186096238E+08 

2=  2443704.500  -. 10649712492385E+09  . 88927994700686E+08 

.39557176239299E+08 

3=  2443718.500  -.98717251872146E+08  .52788739659046E+08 

. 17809766287200E+Q8 


1-2 


4=  2443732.500  -. 79958412568649E+08  . 44646799926066E+08 

. 11065755727868E+08 

5=  2443746.500  - .80061334094818E+08  .62853154790499E+08 

. 21537849807726E+08 

6=  2443760.500  -. 13285945807806E+09  . 73886441590330E+08 

.347 1665 1073760E+08 

1=  time  (Julian  day),  X  (geocentric-equatorial),  Y  (geocentric-equatorial), 
Z  (geocentric-equatorial) 

2=  time  (Julian  day),  X  (geocentric-equatorial),  Y  (geocentric-equatorial), 
Z  (geocentric-equatorial) 

3-6=  same 

The  output  format  for  GEOEQ  is  Format  (1F12.3,  3E30.14)  for  each  record. 
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APPENDIX  J 


PROGRAM  DELAIP 


PROGRAM  DEL ALP l INPUT . OUTPUT , TAPE  1 , TAPE2 , TAPE3 ) 

DIMENSION  T ( 700 ) . ALPHA ( 700 > , DE LT A ( 700 ) , EP0CH(700) 
DIMENSION  TTT ( 700 ) 

REAL  MNA.MND. MENA. MEND 

COMMON/ P LOTD/XMI N , XMAX , YMI N , VMAX . I  CHAR , I  AXIS , PLOTM!  11,51) 
ICT  =  0 

C  READ! 1.91 )EP0CH1 

C  READ12.91 (EP0CH2 

91  FORMAT! 1F12.3) 

DO  50  1*1. 1000 
READ! 1 ,4)TT.X1 . VI .21 
I F ( EOF ( 1 ) ,NE-0.0)GO  TO  60 
READ! 2 , 4 ) TT 1 , X2 , Y2 , Z2 
IF(E0F(2) . NE.OOIQO  TO  60 

4  FORMAT! 1F12.3.3E30. 14) 

IF! ICT.NE.O)GO  TO  51 
I F ( TT . N£ . TT 1 )4 1 , 51 

41  PRINT  42.TT.TT1 

42  FORMAT! »  INITIAL  TIMES  ARE  NOT  EQUIVALENT.  TT  =  *,E30.14, 

1 *TT 1 =  - ,£30. 14) 

GO  TO  90 
51  CONTINUE 

EPOCH! I ) =TT 

T ( I )*EPOCH( I ) -EPOCH! 1 ) 

IF(1. £0.2)15, 14 
15  TDEL*EP0CH(2)-EP0CH( 1) 

14  CONTINUE 

C  ALPH1.ALPH2  ARE  RIGHT  ASCENSIONS 

C  0ELTA1 .DELTA2  ARE  DECLINATIONS 

C  ALPH 1 , ALPH2 .DELTA1 .DELTA2  ARE  RETURNED  IN  ARCMINUTES 

PI *3 . 1415926535898 
ARCMIN- 10800. OO/PI 
ALPH1*ATAN(Y1/X1 ) 

CALL  QUAD! ALPH1 .XI . VI ) 

ALPH 1* ALPH 1*ARCMIN 
ALPH2*AT  AN( Y2/X2 ) 

CALL  OUAD! ALPH2.X2. Y2) 

ALPH2* ALPH2  * ARCMIN 
R01*X1**2+Y1*'2 
RO 1 *SQRT ( RO 1 ) 

R02*X2«*2+Y2-'2 

R02*S0RT(R02) 

DELTA1-ATAN! Z1/R01 ) 

DELTA  1 -DELTA  1* ARCMIN 
DELTA2*AT  AN! Z2/R02 ) 

DELTA2-DELTA2* ARCMIN 
PRINT  18,  ALPH 1 . ALPH2 . DE LT A  1 .DELTA2 
18  FORMAT! 4E30. 14 ) 

C  ALPHA  IS  OIFFERENC  IN  RIGHT  ASCENSIONS  AT  TIME  T 

C  DELTA  IS  DIFFERENC  IN  DECLINATIONS  AT  TIME  T 

ALPHA! I )*ALPH1-ALPH2 
OELTAd )*0ELTA1-DELTA2 
ICT-ICT+1 
TTT( I )*TT1 
50  CONTINUE 

60  CONTINUE 
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I F ( ICT . EQ . 0 ) 10. 11 
PRINT  12 

FORMATC  E-O-F  ENCOUNTERED  AT  BEGINNING  OF  INPUT.  NO  INFORMATION 
1  ON  AT  LEAST  ONE  OF  THE  INPUT  TAPES.*) 

GO  TO  90 
PRINT  7.'ICT 

FORMAT ( ”  ICT*  “.113) 

PRINT  13 


FORMAT (6X, *TIME( JULIAN  DAYS)* . 1 2X . *T { I ) * . 1 3X 
1  . » ALPHA ( ARCMINUTES ) » , 1 3X , *DELT A ( ARCMINUTES ) * ) 
T ( I )  IS  COMPUTED  FOR  PLOTTING  PURPOSES 
MCT  *0 
MNA=0. 0 
MND=0 . 0 
VARA-0 . 0 
VARD=0.0 
DO  40  1*1, ICT 
MNA  =  MNA  +  ALPHA( I  ) 

MND=MND+DELTA( I  ) 

MCT  -MCT  + 1 
T(I)*T(I) /TOE  L 

PRINT  5 , TTT ( I ) . T( I ) . ALPHA ( I ) . DELT A ( I  ) 
WRITE(3.5)TTT( I ) . T ( I ) , ALPHA ( I ) , DELTA ( I ) 

FORMAT ( 1F12.3.3E30. 14) 

CONTINUE 
RMT  -MCT 
MENA-MNA/RMT 
MEND-MND/RMT 
D0 1  dd*1.ICT 

VARA-VARA+ ( ALPHA (dd) -MENA)* *2 
VARD  =  VARD+ ( DELTA ( dd ) -MEND ) *  *2 
CONTINUE 

SIGA=SORT(VARA/RMT) 

SIGD*SORT( VARO/RMT) 

PRINT  8 . MENA . MEND , SIGA . SIGD 
AMINV*AMINE(T.ICT) 

AMAXV*AMAXE(T. ICT) 

XMIN*AMINV 

XMAX*AMAXV 

AMINV- AMINE (ALPHA, ICT) 

AMAXV-AMAXE ( ALPHA , ICT ) 

YMIN*AMINV 
YMAX* AMAXV 
PRINT  9.YMIN.YMAX 
FORMAT ( 1H.4E30. 14 ) 

AMINV* AMINE (DELTA , ICT ) 

AMAXV*AMAXE( DELTA , ICT ) 

YMIN-AMINV 
YMAX*AMAXV 
PRINT  9.YMIN.YMAX 
FORMAT ( 1H.2E30. 14) 

CONTINUE 

STOP 

END 

SUBROUTINE  QUAD( ALPHA, X.Y) 

PI *3 . 1415926535898 
IF(ALPHA.GT. 0.0)1 ,2 
IF(Y.GT.O.O)  ALPHA* ALPHA 
IF(Y.LT.O.O)  ALPHA “ALPHA+P I 
RETURN 

IF(ALPHA.LT.O.O)  3.4 
IF(Y.GT.O.O)  ALPHA-ALPHA+PI 
IF(Y.LT.O.O)  ALPHA*ALPHA+2 .0*PI 
RETURN 

IF(ALPHA.EO.O.O)  5.6 
IF(X.GT.O.O)  ALPHA-0.0 
IF(X.LT.O.O)  ALPHA-PI 
CONTINUE 
RETURN 
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PROGRAM  DELAPM( INPUT . OUTPUT , TAPE  1 . TAPE2 , TAPE3  ) 

DIMENSION  T(700) .ALPHA (700) .DELTA (700) ,EP0CH(700) 

DIMENSION  TTT ( 700 ) 

REAL  MNA . MND , MENA . MENO 

COMMON/PLOTD/XMIN , XMAX . YMIN , YMAX , ICHAR . I  AXIS . PLOTM( 1 1 .61 ) 

ICT-0 

C  READ( 1.91 ) EPOCH  1 

C  R£AD( 2.91 ) EP0CH2 

91  FORMAT ( IF  12 . 3 ) 

00  50  I«1, 1000 
READ( 1 ,4)TT ,X1 , Y1 .21 
I F ( EOF ( 1 ) ,NE.O.O)GO  TO  60 
RE AO ( 2 , 6 ) ALPH2 , DELT  A2 
IF(EOF(2).NE.O.O)GO  TO  60 

4  FORMAT ( 1 F 1 2 . 3 . 3E30 .14) 

6  F0RMAT(2E30. 14 ) 

EPOCH( I ) -TT 

T ( I  )*EPOCH(I )-EPOCH( 1 ) 

IF(I.E0.2)15, 14 
15  TDEL-EPQCH( 2 ) -  EPOCH (  1 ) 

14  CONTINUE 

C  ALPH1.ALPH2  ARE  RIGHT  ASCENSIONS 

C  DELTA  1 , 0ELTA2  ARE  DECLINATIONS 

C  ALPH1 . ALPH2.DELTA1 .DELTA2  ARE  RETURNED  IN  ARCMINUTES 

PI *3 .  1415926535898 
ARCMIN- 10800.00/PI 
ALPH1 »ATAN( Y 1/X1 ) 

CALL  OUAD( ALPH 1 . X 1 , Y 1 ) 

ALPH1 -ALPH1 -ARCMIN 
R01-X1--2+Y1--2 
R01«SQRT(R01 ) 

DELTA1=ATAN(Z1/R01 I 
DELTA 1-DELTA1* ARCMIN 
18  FORMAT ( 4E30 . 14 ) 

C  ALPHA  IS  DIFFERENCE  IN  RIGHT  ASCENSIONS  AT  TIME  T 

C  DELTA  IS  DIFFERENCE  IN  OECLINATIONS  AT  TIME  T 

ALPHA( I J-ALPH1 -ALPH2 
DELTA ( I  )»DELTA1-DELTA2 
ICT-ICT+1 
TTT ( I )=TT 
50  CONTINUE 

60  CONTINUE 

I F ( I CT . EO .0)10. 11 

10  PRINT  12 

12  FORMAT ( *  E-O-F  ENCOUNTERED  AT  BEGINNING  OF  INPUT.  NO  INFORMATION 
1  ON  AT  LEAST  ONE  OF  THE  INPUT  TAPES.-) 

GO  TO  90 

1 1  PRINT  7 , ICT 

7  FORMAT ("  ICT-  " ,113) 

C  T ( X )  IS  COMPUTEO  FOR  PLOTTING  PURPOSES 

MCT-0 
MNA-O. 0 
MND -0.0 
VARA-0 -  0 
VARD-0.0 
DO  40  1*1  . ICT 
MNA-MNA*ALPHA( I ) 

MND-MNO+DELTA( I ) 

MCT-MCT+ 1 

T( I  )»T( 1  J/TDEL 

WRITEO,  5)TTT(  I  )  ,T(I  )  .ALPHA  (  I).DELTA(I) 

5  FORMAT ( 1 F 1 2 . 3 . 3E30 . 14 ) 

40  CONTINUE 

RMT-MCT 
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MENA*MNA/RMT 
MEND*MND/RMT 
D0 1  JU= 1 . ICT 

VARA=VARA+( ALPHA (UJ)-MENA)**2 
VARO  =  VARD+( DELTA! UJ ) -MEND  )-*2 
1  CONTINUE 

SI  GA  =  SORT ( VARA/RMT ) 

SI GD'SQRT ( VARD/RMT  ) 

PRINT  E. MENA. MEND, S1GA.SIGD 
AM I NV- AMINE ( T . ICT ) 

AMAXV*AMAXE(T . ICT ) 

XMIN*AMINV 

XMAX-AMAXV 

AMINV* AMINE ( ALPHA , ICT) 
AMAXV=AMAXE ( ALPHA . ICT) 

YMIN* AMINV 
YMAX* AMAXV 
PRINT  9.YMIN.YMAX 

8  FORMAT! 1H.4E30. 14) 
AMINV*AMINE!D£LTA, ICT) 
AMAXV«AMAXE( DELTA . ICT ) 
YMIN=AMINV 

YMAX= AMAXV 
PRINT  9.YMIN.YMAX 

9  FORMAT! 1H, 2E30. 14 ) 

90  CONTINUE 

STOP 

END 

SUBROUTINE  OUAD 1  ALPHA . X , Y  ) 

PI *3 . 1415926535898 
IF! ALPHA. GT. 0.0)1. 2 

1  I F ( Y . GT .0.0)  ALPHA x ALPHA 
IF(Y.LT.O.O)  ALPHA'ALPHA+PI 
RETURN 

2  IF( ALPHA . LT .0.0)  3.4 

3  IF(Y.GT.O.O)  ALPHA-ALPHA+PI 
IF(Y.LT.O.O)  ALPHA«=ALPHA  +  2 .0*PI 
RETURN 

4  IF! ALPHA. EO. 0.0)  5.6 

5  IF(X.GT.O.O)  ALPHA=0 . 0 
IF(X.LT.O.O)  ALPHA*PI 

6  CONTINUE 
RETURN 
END 


J-4 


NSWC  TR  83-43 


SAMPLE  INPUT  FOR  TAPE  1 


These  are  JPL  coordinates  for  planet  "Mercury.” 


Time 


1=  2443690.500  - . 77964575183693E+08 

.68152401023670E+08 

2=  2443704.500  - . 10649895546677E+09 

. 39561718564281E+08 

3=  2443718.500  - .98717650891779E+08 

. 17814249153166E+08 

4=  2443732.500  - . 79964288575854E+08 

. 11063858710573E+08 

5=  2443746.500  - . 80059614384767E+08 

•21535914238598E+08 

6=  2443760.500  - . 13286136735694E+09 

.34712225207949E+08 


. 14285139938627E+09 
. 88922413933842E+08 
.52787733471752E+08 
. 44635600380581E+08 
. 6284972545 1473E+08 
. 73876199041347E+08 


1=  time  (Julian  days),  X  (geocentric-equatorial),  Y  (geocentric-equatorial), 
Z  (geocentric-equatorial) 


2-6=  same 


SAMPLE  INPUT  FOR  TAPE  2 

These  are  Kepler  generated  coordinates  for  planet  "Mercury." 


Time 


1=  2443690.500  - . 77974525023110E+08 

.68147186096238E+08 

2=  2443704.500  - . 10649712492385E+09 

. 39557 176239299E+08 

3=  2443718.500  - .98717251872146E+08 

. 17809766287200E+08 

4=  2443732.500  - . 79958412568649E+08 

. 11065755727868E+08 

5=  2443746.500  - . 80061334094818E+08 

. 21537849807726E+08 

6=  2443760.500  - . 13285956807806E+09 

.347 1665 107376E+08 


.  14285226329953E+09 
.  88927994700686E+08 
.52788739659046E+08 
.  44646799926066E+08 
.  62853154790499E+08 
.  73886441590330E+08 


1=  time  (Julian  days),  X  (geocentric-equatorial),  Y  (geocentric-equatorial), 
Z  (geocentric-equatorial) 


2-6=  same 


The  input  format  for  the  original  program  DELALP  is  Format  (1F12.3,  3E30.14)  as  in 
records  one  through  six  for  both  JPL  and  Kepler  generated  coordinates,  Tape  1  and 
2.  The  input  formats  for  the  revised  program  DELAPM  are  Format  (1F12.3,  3E30.14) 
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»r  the  time  and  position  histories  generated  by  JPL,  Tape  1  and  Format  (2E30.14) 
»r  the  right  ascension  and  declination  angles  generated  by  program  select,  Tape 


IMPLE  OUTPUT  FOR  PROGRAM  DELALP 

Statistical  output  generated  by  the  print  statement  for  planet  "Mercury." 
:T=  562 

. 47901140479242E-03  .91930569345719E-03  .17183945 

L7541E+00  . 7 1253067327 114E-01 

- . 71675558801508E+00  . 73359012082801E+00 

- . 2525446532499 7E+00  . 233431 14920863E+00 

ICT  =  562  ■*  number  of  position  histories 

mean  -  right  ascension,  mean  declination,  standard  deviation  - 
ight  ascension,  standard  deviation  -  declination 

minimum  value  -  right  ascension,  maximum  value  -  right  ascension 
minimum  value  -  declination,  maximum  value  -  declination 

Sample  output  for  Tape  3  which  reflects  the  changes  in  angles. 

=  2443690.500  0. 

. 13530473608989E+00 

=  2443704.500  . 10000000000000E+01 

. 11827692166844E+00 

=  2443718.500  .20000000000000E+01 

. 13485020775261E+00 

=  2443732.500  . 30000000000000E+01 

-.68717342963282E-01 

=  2443746.500  . 40000000000000E+01 

- . 38834355153085E-01 

=  2443760.500  .50000000000000E+01 

- . 78883304417104E-01 

1=  time  (Julian  days),  counter,  A  right  ascension 
A  declination 

2=  time  (Julian  days),  counter,  A  right  ascension 
A  declination 

3-6=  same 

he  output  formats  for  the  printed  statistics  are  Format  ("ICT=",  113)  for  the 
umber  of  position  histories  in  the  first  line.  Format  (1H,  4E30.14)  for  mean 
nd  standard  deviation  in  lines  2  and  3  and  Format  (1H,  2E30.14)  for  minimum 
nd  maximum  values  in  lines  4  and  5.  The  output  format  for  the  tape  is  Format 
1F12.3,  3E30.14).  This  format  is  repeated  for  each  record. 


- . 17574213532498E+00 
. 13521342672175E+00 
. 33026276330929E-01 
.47460765269352E+00 
.552390582 16 102E-01 
. 22341 148892883E+00 


PROGRAM  SELECT ( INPUT .OUTPUT .TAPE  1 . TAPE2 . TAPE3 . T APE4 . TAPE5 , 

6TAPE6 . T  APE  1 1  ) 

THESE  FACTORS  CONVERT  THE  GIVEN  VALUES  FROM 
FACT  -  ARC  SECONDS  TO  RADIANS 
FACTi  -  TIME  SECONDS  TO  RADIANS 
F ACT 2  -  DEGREES  TO  RADIANS 

F ACT 3  -  TIME  HOURS  TO  RADIANS 

REAL  IN.LN 

DIMENSION  RAl  8 ) ,DEC( 8 ) ,DXP( 8 ) ,DVP( 8 ) ,D2P( 8 ) . RMAG!  8 ) 

DIMENSION  AD (6 ) . EC (6 ) . INI  6 ) . LN(6 ) . OMEGA! 6 ) . AMO! 6 ) . AM  1 ( 6 ) . AM! 6 ) 

DATA  FACT/  484813681 10953E-05/ , FACTi/ . 72722052 166429E -04/ 

DATA  FACT2/ . 1 74532925 1 9943E -01 /. FACT3/ . 26 1 799387799 1 4/ 

DATA  AD/ 1 . 000000003. . 3870986732 . .7233322265 . 1 . 52369076,5.203070156 
A  .9.538509603/ 

DATA  EC/. 01671532949782. .20562942641 17, .006775929468209, 

B  .09339472415910. .04825983079881 . .05315110813689/ 

DATA  I N/ . 0000898955 1 986890 . . 1222021384653, .05923236598127, 

C  .03223176504070. .02280397821074. .04342907360759/ 

DATA  OMEGA/ 1 . 783756851 , 1 . 339368463 , 2 . 28392598 , 5 . 852403589 . 

D  .2583043378.1.599695040/ 

DATA  LN/3. 049977 190305. .8323191351502, 1.328556712384, 

E  .8562002543165. 1 .746015566160. 1 .974220936589/ 

DATA  AMO/3 . 073477677683 . 1 . 248550488749 , 1 . 160433857370. 

F  3.899288800052. 1 .505867745775, .9566213166208/ 

DATA  AMI/. 0172021247, .0714247869, .027962452, .0091461066. 

G  .001450104, .0005840133/ 

DATA  TEP0CH/244369O. 5/ . DELT/ 14.0/ 

DATA  TEND/2451558. 5/ 

DATA  PI/3.1415926535898/ 

DATA  TU/58. 13244087/ 

TIME*TEPOCH 
EDAY-TIME-TEPOCH 
DO  30  K* 1 ,6 

AM(K)«AMO!K)+AM1(K)*EDAY 

continue 

DO  40  KZ»1 .6 

CALL  E  LC0RD !KZ,A0!KZ),EC!KZ), IN! KZ ) . AM( KZ } , OMEGA ( KZ ) . LN< KZ ). X , Y , Z ) 

IFtKZ.EO. 1 )  GO  TO  20 

OXP ( KZ ) *OX*X -XE 

DYP! KZ )-DY*Y-YE 

DZP ( KZ  l-OZ-Z-ZE 

ARCMIN- 10800.00/PI 

RMAG ( KZ ) *  SORT ( DX-DX+D Y-DY+DZ-DZ ) 
RA(KZ)«ARTNO(DY ,DX )» ARCMIN 
DEC(KZ)*ASIN( DZ/RMAG! KZ ) ) -ARCMIN 
WR I TE ( KZ , 7202 )RA(KZ) .DEC(KZ) 

RA(KZ)*ARTN0(DY,DX)/FACT3 
DECIKZ  >  =  AS  IN! DZ/RMAG! KZ ) J/FACT2 
GO  TO  40 
20  DXP! 1 >*XE*X 
DYP!  1  ) *  Y E  *  Y 
DZP!  1 )*ZE*Z 
DXP ( 8  ) * -XE 
OYP ( 8 ) ■ - Y  £ 

DZP! 8 )*-ZE 

RMAG! 8 ) -SORT (XE-XE+YE-YE+ZE-ZE  ) 

RA ( 8 ) * ARTNO! -YE . -XE I/FACT3 
DEC! 8 )*ASIN( - ZE/RMAG! 8) )/FACT2 


n  o  n  n 


r 
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40  CONTINUE 

AM  350.737486+12 .  190749 19 14 -DA VS ) 'FACT2 
B- ( 296 . 104608+1 3 . 0649924465 -DAYS )'FACT2 
C*( 358. 475833+0. 9856002669-DA VS )'FACT2 
DM 1 1 . 250889+13. 229350449'DAYS ) 'FACT2 
T  A  3  2 . *A 

CB=COS(B)  i  CB2A=C0S(B-TA)  $  C2A=C0S(TA> 

SBrSIN( B )  $  SB2 A  =  S IN( B - T A )  $  S2A=SIN(TA)  $S2B*SIN( 2 . »B ) 

SC-SIN(C)  $  SD-SIN(D)  $  SBPD-SIN( B+D )  $  SD2A-SIN( D-TA ) 

SDB=SIN(0-B ) 

RL  -RMLM+ (22639. 5+ SB -4586 . 426+SB2A+2369 . 902+S2A+769 . 016-S2B-668 . 1 
1  1 1 -SC  I 'FACT 

RP  M 18461 .48-SD+1010. 1 8 • SBPD -999 . 695'SDB -623 . 658 *SD2 A ) * F ACT 
PARMM  3422 . 7+ 186 . 5398*CB+34 . 3 1 1 7 +  CB2 A+28 . 2333-C2A ) 'FACT 
CRL-COS(KL)  $  SRL-SIN(RL) 

CRP-COS(RP)  %  SRP*SIN( RP )  $  SPARM-SINI PARM ) 
AI-CRP'SRL'CEPS-SRP'SEPS  $  A2-CRP'CRL 
ALPHAM-ARTNQI  A  1 , A2  ) 

SDELTM»CRP'SRL*SEPS+SRP*CEPS 
CDELTM-SORT( 1 . -SDELTM'SOELTM) 

HRANG-TS-ALPHAM 

X-CDELTM'SIN(HRANG) 

Y-CDELTM*C0S(HRANG)-CLA*SPARM 
Z-SDELTM-SLA-SPARM 
RH0-SQRT(X-X+Y»Y+Z-Z ) 

S0 1 -Z/RHO  $  CD1=S0RT( 1 . -SDI'SDl ) 

E*ARTNQ( X . Y )  $  TSE-TS-E 

RA(7)-TSE/FACT3  $  DEC(7)-ASIN(SD1 )/FACT2 
I F ( RA ( 7 ) . LT . O . )  RA(7)«RA(7)+24. 

IF  ( RA ( 7 ) . GE . 24 .  )  RA ( 7 ) -RA ( 7 ) - 24 . 

DXP ( 7 ) -CD  1 -COS ( TSE ) 

DYP(7)-CD1*SIN(TSE) 

OZP (  7  I-S0 1 

RMAG( 7 ) -SORT ( DXP ( 7 ) -DXP (7)+DYP(7)*DYP(7 )+DZP( 7 )-DZP(7) ) 

PRINT  50 . SP ARM , TS . E . X , Y . Z 
PRINT  50.TA.B.C.D.RL.RP 

PRINT  50 . CRL , SRL . CRP . SRP , ALPHAM . A  1 . A2 , TSE . SD 1 , CD1 
50  FORMAT ( 6E 19 , 12) 

7202  FORMAT(2E30. 14) 

Y-Y+CLA'SPARM  $  Z-Z+SLA-SPARM 
PRINT  50.X.Y.Z 
TIME-TIME+DELT 
I F ( T IME . LT . TEND )G0  TO  16 
END 

SUBROUTINE  ELCORD ( KPL . A . E , RI . RL , RO , RH , X , Y . Z ) 

METHOD  IN  SMART  USING  LONGITUDE  OF  PERIHELION  AND  TRUE  ANOMALY 
INSTEAD  OF  ARGUMENT  OF  PERIHELION  AS  IN  DR  0[S  PROGRAM 

DIMENSION  RNEW( 3 ) , ROLD( 3 ) 

DATA  SORO/- 1 . 7202 10850E -02/ , TOL/ 1 .745E-06/ 

DATA  EPS/ .4092061859/ 

EO*RL+E*SIN(RL )+.5'E-E-SIN(2.-RL) 

15  DELEOM  RL-E0+E*SIN(E0) )/( 1 . -E*COS(  EO) ) 

EO-EO+DELEO 

IF  (ABS(DELEO).LT.TOL)  GO  TO  25 
GO  TO  15 
25  TE-TAN( EO' . 5 ) 

V-2  *ARTNQ(TE-SQRT( 1 ,+E  )  . SORT ( 1 . -E ) ) 

R-A-( 1 . -E*COS( EO) ) 

SH-SIN(RH)  $  CH-COS(RH)  $  ANG-RO-RH+V  $  SA-SIN(ANG)  $  CA-COS(ANG) 
CI-COS(RI)  $  A11-CH-CA -SH'SA-CI  $  A2  1 -SH'CA+CH* SA «CI 
ROLD (  1  ) -R' A  1 1  $  ROLO( 2  )*R'A2 1  *  R0lD( 3 )-R*SA»SIN( RI ) 

RNEWl 1 ) -ROLD ( 1 ) 

RNEW12  )*ROLD( 2 ) -COS (EPS )-SIN( EPS)*R0LD(3) 

RNEWf 3  )*SIN( EPS  )'ROLD( 2 )+COS( EPS ) -ROLDf  3 ) 

X-RNEW(I)  $  Y-RNEW(2)  $  Z-RNEW13) 

RETURN 
END 


K-2 
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This  subroutine  was  extracted  from  a  F/C  program.  In  our  exercise,  not  all 
iables  referenced  in  this  routine  were  utilized.  Therefore  several  variables 
r  appear  undefined. 

The  input  items  for  program  SELECT  are  contained  in  data  statements  within 
:  program. 


1PLE  OUTPUT  FOR  TAPE  2  PLANET  MERCURY 
Right  Ascension 

=  . 71 176495109891E+04 
=  . 84082320444799E+04 
=  . 91 1186979341 13E+04 
=  -90493253525415E+04 
=  . 85 1 19521362692E+04 
=  . 90552286856794E+04 
=  . 10445917609449E+05 
=  . 11797161582502E+05 
=  .13048641098172E+05 
=  . 14272139430453E+05 
=  . 15350209710689E+05 


Declination 

.  13632289664085E+04 
.95480166004741E+03 
.  54237762899285E+03 
.  41338946820392E+03 
■71685273096590E+03 
.77182633642678E+03 
.  26981732572514E+03 
.  37956019556662E+03 
.  94346874995080E+03 
.  13475293043596E+04 
.  15233571132286E+04 


1=  right  ascension-angle  (arc  minutes),  declination-angle  (arc  minutes) 

11=  right  ascension-angle  (arc  minutes),  declination-angle  (arc  minutes) 

-  after  a  change  in  time  -  -  after  a  change  in  time  - 

e  output  format  for  program  SELECT  is  Format  (2E30.14).  This  format  is  repeated 
r  each  record. 
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APPENDIX  L 
PROGRAM  PLOT 


PROGRAM  PLOT!  INPUT .OUTPUT . TAPE  1  . 

DIMENSION  T ( 600  I  ,C(600)  . ALPHA(600) 

DIMENSION  DELTA! 700 ) .DUM! 16).TITlE(3I.EXXLAB(2) 
DATA  EXXLAB/ 1 1HARC  MINUTES/ 

DATA  EXlABB/ 10HFORTNIGHTS/ 

READ  15. TITLE 
READ  -,KC 
READ  *,L 
READ  -.N 
DC  10  1=1. 562 . 5 
J=( 1*4  )/5 

READ! 1,5)T(U).C(J) , ALPHA! J ) .DELTA! J ) 

IF! EOF( 11)20.8 
8  READ  (1,5)  DUM 

K  =  U 

10  CONTINUE 

20  CALL  XRANGEIO. 0.600.0) 

CALL  VRANGE! -2. 0,2. 5) 

IFIN.EO. 1 130.40 
30  CALL  PLOTT!C. ALPHA. K) 

GO  TO  50 

40  CALL  PLOTT(C. DELTA. Kl 

50  CALL  LABELL! 30, L ,KC . TITLE ) 

CALL  LABELL! 20. 5 .- 1 1 . EXXLAB ) 

CALL  LABELH2. 31. 10. EXLABB) 

CALL  PAUSEE(ICH) 

5  FORMAT! iF 12 . 3, 3E30. 14 ) 

15  FORMAT ( 3A 10 ) 

END 


SAMPLE  INPUT  FOR  TAPE  1  PLANET  MERCURY 

This  input  represents  the  differences 
and  declination  angles. 

1=  2443690.500  0. 

. 13530473608989E+00 

2=  2443704.500  . 10000000000000E+01 

. ii827692l66844E+00 

3=  2443713.500  . 20000000000000E+01 

. 13485020775261E+00 

4=  2443732.500  . 30000000000000E+01 

-. 68717342963282E-01 

5=  2443746.500  .40000000000000E+01 

- . 38834355153085E-01 

6=  2443760.500  . 50000000000000E+01 

-. 78883304417104E-01 


between  Kepler  and  JPL  right  ascension 

-.17574213532498E+00 
.13521342672175E+00 
. 33026276330929E-01 
.47460765269352E+00 
.55239058216102E-01 
22341 148892883E+00 


1=  time  (Julian  days),  Counter  t(I),  A  right  ascension 
A  declination 
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2=  time  (Julian  days),  Counter  t(I),  A  right  ascension 
A  declination 

3-6  =  same 

The  input  format  for  the  tape  is  Format  (1F12.3,  3E30.14).  This  format  is  re¬ 
peated  for  each  record. 

The  following  is  sample  input  for  program  PLOT  to  be  inserted  during  the 
time  of  execution  at  the  terminal. 

Read  Input 


Title 

KC 

L 

N 


Mercury  -  Right  Ascension 
25 
24 
1 


Title  is  the  title  of  the  plot  and  is  formatted  Format  (3A10).  KC  is  the  num¬ 
ber  of  characters  in  the  title.  L  is  the  column  in  which  the  title  should  start 
to  be  written.  N  is  the  control  parameter  that  decides  if  right  ascension  or  de¬ 
clination  is  to  be  plotted.  If  N  equals  one,  right  ascension  is  plotted.  For  any 
other  number,  declination  is  plotted.  KC,  L,  and  N  are  all  free-formatted  input 
items . 
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